DIME/dimemcv1.05.f

   1       subroutine dimeinit
   2
   3       implicit double precision(a-y)
   4       implicit complex*16(z)
   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
  10      &,hh,lmax
  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
  13       integer iin,nch
  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
  16       common/ipars/nch
  17       common/mom/q
  18       common/int/ip
  19       common/mompt/pt1,pt2,ptx
  20       common/vars/s,rts,mmes,yx
  21 ccc new
  22       common/vars0/ mf127, mf1525, m0, mmes0, mmes1, mmes2, mp,
  23      &              mwidth, pi, rt2, ebeam, sum, sum1, weightm,
  24      &              bjac, bp
  25       common /ivars/ ncut, ll, icut, nev, num
  26 ccc new
  27       common/cuts/etaelmax,etaelmin,ptelmin,ptphmin,ecut,rmax,rmin,mcut
  28       common/flags/iin, pflag, fsi, ppbar, output, cuts, unw
  29       common/gluons/g1,g2,kt1,kt2
  30       common/levicivita/epsilon
  31       common/ang/cost,costa,costb
  32       common/alph/alfas
  33       common/hvars/sh,th,uh
  34       common/wvars/sig0,bb,t1,t2
  35       common/alphas/alphap,alpha0,alphapr,alpha0r,alphapm,alpha0m
  36       common/sec/cpom,cf,crho,aff,ar
  37       integer iii
  38       common/ivar/iii
  39       common/ff/formf
  40       common/ffpars/bexp,ao,bo,aoo,boo
  41       common/regge/mregge
  42 ccccc hepevt output
  43       integer nmxhep,kk
  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
  50       INTEGER MAXNUP
  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),
  57      &              SPINUP(MAXNUP)
  58 ccccc rambo variables
  59       integer npart,nrun
  60       double precision pin(4),am(100),pout(4,100),wt,ein
  61       common/ramboc/ pin, am, pout, wt, ein,
  62      &npart, nrun
  63 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  64 c                                                         c
  65 c     Dime MC for the central exclusive production of     c
  66 c     meson pairs by double Pomeron exchange:             c
  67 c                                                         c
  68 c     p(1) p(2) --> p(3) + M_1(5) M_2(6) + p(4)           c
  69 c                                                         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
  72 c     component, with:                                    c
  73 c                                                         c
  74 c     1,2 = transverse components                         c
  75 c     3   = beam component                                c
  76 c     4   = energy                                        c
  77 c                                                         c
  78 c     PDG number for ith particle in arrary ID(i)         c
  79 c                                                         c
  80 c     Also gives event record in HEPEVT or LHE format     c
  81 c     (others are available upon request)                 c
  82 c                                                         c
  83 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  84 c                                                         c
  85 c     Particles generated:                                c
  86 c                                                         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
  92 c                                                         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
  96 c     width included                                      c
  97 c                                                         c
  98 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  99 c                                                         c
 100 c     User defined cuts can readily be implemented        c
 101 c     in subroutine 'icut'                                c
 102 c                                                         c
 103 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 104 c                                                         c
 105 c     This is version 1.05           : 17 Sep 2014        c
 106 c                                                         c
 107 c     Comments etc to: lucian.harland-lang@durham.ac.uk   c
 108 c                                                         c
 109 c     If you use the code please reference (and see for   c
 110 c     details of model used):                             c
 111 c                                                         c
 112 c     L.A. Harland-Lang, V.A. Khoze, M.G. Ryskin          c
 113 c     and W.J. Stirling arXiv:1105.1626                   c
 114 c                                                         c
 115 c     L.A. Harland-Lang, V.A. Khoze, M.G. Ryskin          c
 116 c     and W.J. Stirling arXiv:1204.4803                   c
 117 c                                                         c
 118 c     L.A. Harland-Lang, V.A. Khoze and  M.G. Ryskin      c
 119 c     arXiv:1312.4553                                     c
 120 c                                                         c
 121 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 122
 123 ccAM      open(35,file='exerec1.dat') ! event record
 124
 125 ccAM      rts=1.96d3         ! centre of mass energy
 126       write(6,*) "energy", rts
 127 cccc  Some basic cuts, imposed in subtroutine 'icut'. Other user defined cuts can readily be implemented in subroutine
 128 cccc  note: in rhorho case these cuts are imposed on rho's, and not their decay productions. Cuts on the decay products
 129 cccc  can be imposed by editing 'icut'
 130
 131 ccAM      rmax=1.8d1        ! max meson rapidity
 132 ccAM      rmin=-1.8d1       ! min meson rapidity
 133 ccAM      ecut=2d0          ! min meson p_t
 134
 135 c -- physics parameters
 136
 137 ccAM      pflag='rho'    ! Process generated - see preamble for options
 138 ccAM      fsi='true'      ! phenomenological model for exclusive suppression in Pom Pom --> meson pair. To turn on/off --> 'true/false'
 139 ccAM      formf='orexp'   ! meson - pomeron form factor.
 140 ccAM      ppbar='false'   ! set true if ppbar collisions
 141           output='lhe' ! Event record style, HEPEVT or LHE
 142 ccAM      cuts='true'     ! Impose cuts or not
 143 ccAM      unw='true'      ! Set = 'true' for unweighted events
 144 ccAM      iin=1           ! Model for soft survival factor, as described in arXiv:1306.2149. Default = 1
 145
 146 cccccccc
 147
 148 ccAM      ntotal=1000000            ! no. of runs for weighted events
 149 ccAM      nev=100                  ! no. of unweighted events generated to event record
 150
 151 ccccccccc
 152
 153       if(formf.eq.'exp')then           ! Set parameters for Pomeron-Meson form factor in function 'fpi(x)'
 154          bexp=1d0/2.2d0
 155       elseif(formf.eq.'orexp')then
 156          bo=1d0/1.1d0
 157          ao=dsqrt(0.5d0)
 158       elseif(formf.eq.'power')then
 159          aoo=1.7d0
 160       endif
 161
 162 ccccccccccccccccccccccccccccccccccccccc
 163
 164 ccccccccc   Start of main body of code
 165
 166 cccccccccccccccccccccccccccccccccccccc
 167
 168       call initpars(iin)   ! Initialise soft survival parameters
 169       call calcop          ! proton opacity
 170       call calcscreen      ! screening amplitude
 171
 172       if(pflag.eq.'pipm')then
 173          mmes=0.13957018d0      ! pi+/- mass, PDG 2011 value
 174          sig0=13.63d0
 175          alphapm=0.7d0
 176          alpha0m=-0.7d0*mmes**2
 177          cf=31.79d0/0.389d0
 178          crho=4.23d0/0.389d0
 179
 180          nup=6
 181          istup(5)=1
 182          idup(5)=211
 183          mothup(1,5)=1
 184          mothup(2,5)=2
 185          icolup(1,5)=0
 186          icolup(2,5)=0
 187          vtimup(5)=0
 188          spinup(5)=9.
 189
 190          istup(6)=1
 191          idup(6)=-211
 192          mothup(1,6)=1
 193          mothup(2,6)=2
 194          icolup(1,6)=0
 195          icolup(2,6)=0
 196          vtimup(6)=0
 197          spinup(6)=9.
 198
 199       elseif(pflag.eq.'pi0')then
 200          mmes=0.1349766d0       ! pi0 mass, PDG 2011 value
 201          sig0=13.63d0
 202          alphapm=0.7d0
 203          alpha0m=-0.7d0*mmes**2
 204          cf=31.79d0/0.389d0
 205          crho=0d0
 206
 207          nup=6
 208          istup(5)=1
 209          idup(5)=111
 210          mothup(1,5)=1
 211          mothup(2,5)=2
 212          icolup(1,5)=0
 213          icolup(2,5)=0
 214          vtimup(5)=0
 215          spinup(5)=9.
 216
 217          istup(6)=1
 218          idup(6)=-111
 219          mothup(1,6)=1
 220          mothup(2,6)=2
 221          icolup(1,6)=0
 222          icolup(2,6)=0
 223          vtimup(6)=0
 224          spinup(6)=9.
 225
 226       elseif(pflag.eq.'kpkm')then
 227          mmes=0.493677d0         ! K+/- mass, PDG 2011 value
 228          sig0=11.82d0
 229          cf=17.255d0/0.389d0
 230          crho=9.105d0/0.389d0
 231
 232          nup=6
 233          istup(5)=1
 234          idup(5)=321
 235          mothup(1,5)=1
 236          mothup(2,5)=2
 237          icolup(1,5)=0
 238          icolup(2,5)=0
 239          vtimup(5)=0
 240          spinup(5)=9.
 241
 242          istup(6)=1
 243          idup(6)=-321
 244          mothup(1,6)=1
 245          mothup(2,6)=2
 246          icolup(1,6)=0
 247          icolup(2,6)=0
 248          vtimup(6)=0
 249          spinup(6)=9.
 250
 251       elseif(pflag.eq.'ks')then
 252          mmes=0.497614d0         ! K_0 mass, PDG 2011 value
 253          sig0=11.82d0
 254          cf=17.255d0/0.389d0
 255          crho=0d0
 256
 257          nup=6
 258          istup(5)=1
 259          idup(5)=310
 260          mothup(1,5)=1
 261          mothup(2,5)=2
 262          icolup(1,5)=0
 263          icolup(2,5)=0
 264          vtimup(5)=0
 265          spinup(5)=9.
 266
 267          istup(6)=1
 268          idup(6)=310
 269          mothup(1,6)=1
 270          mothup(2,6)=2
 271          icolup(1,6)=0
 272          icolup(2,6)=0
 273          vtimup(6)=0
 274          spinup(6)=9.
 275
 276       elseif(pflag.eq.'rho')then
 277          mmes0=0.77549d0          ! rho mass, PDG 2013 value
 278          mwidth=0.1491d0         ! rho width PDG 2013 value
 279          sig0=10d0
 280          cf=0d0
 281          crho=0d0
 282
 283          nup=10
 284          istup(5)=2
 285          idup(5)=113
 286          mothup(1,5)=1
 287          mothup(2,5)=2
 288          icolup(1,5)=0
 289          icolup(2,5)=0
 290          vtimup(5)=0
 291          spinup(5)=9.
 292
 293          istup(6)=2
 294          idup(6)=113
 295          mothup(1,6)=1
 296          mothup(2,6)=2
 297          icolup(1,6)=0
 298          icolup(2,6)=0
 299          vtimup(6)=0
 300          spinup(6)=9
 301
 302          do k=7,10
 303             istup(k)=1
 304             mothup(2,k)=0
 305             icolup(1,k)=0
 306             icolup(2,k)=0
 307             vtimup(k)=0
 308             spinup(k)=9.
 309          enddo
 310
 311          idup(7)=211
 312          idup(8)=-211
 313          idup(9)=211
 314          idup(10)=-211
 315          mothup(1,7)=5
 316          mothup(1,8)=5
 317          mothup(1,9)=6
 318          mothup(1,10)=6
 319
 320       endif
 321
 322 c -- other parameters
 323
 324       ebeam=rts/2d0
 325       s=rts*rts
 326       zi=(0d0,1d0)
 327       rt2=dsqrt(2d0)
 328       pi=dacos(-1d0)
 329       bp=rts/dsqrt(2.0d0)
 330       mp=0.93827d0
 331       beta=dsqrt(1d0-4d0*mp**2/s)
 332       s0=1d0
 333
 334 c     Pomeron + t-slope
 335
 336       bb=4d0
 337       bjac=6d0
 338       bjac1=2d0
 339
 340       alphap=0.25d0    ! D-L '92 fit
 341       alpha0=1.0808d0
 342       alphapr=0.93d0
 343       alpha0r=0.55d0
 344
 345       mf127=1.275d0
 346       mf1525=1.525d0
 347
 348       cpom=sig0/0.389d0
 349       aff=-0.860895d0
 350       ar=-1.16158d0
 351
 352       mmes1=mmes
 353       mmes2=mmes
 354
 355 ccccc Initialise rambo (rho0 decay)
 356
 357       if(pflag.eq.'rho')then
 358
 359          mmes=mmes0
 360          npart=2
 361          do j=1,4
 362             am(j)=0.13957018d0
 363          enddo
 364       endif
 365
 366 ccccccc
 367
 368 c     initialise counters etc
 369
 370 c      do hh=1,20
 371
 372       nhist=1
 373       sum=0d0
 374       sum1=0d0
 375       ncut=0
 376
 377       weightm=0d0
 378
 379 c     initialise histograms
 380
 381       do i=1,nhist
 382          call histo3(i)
 383       enddo
 384
 385       num=0
 386
 387 cccccccccccccccccccccccccccccccccccccccccccccccccccccc
 388
 389 c     Incoming protons to event record array
 390
 391 cccccccccccccccccccccccccccccccccccccccccccccccccccccc
 392
 393       ID(1)=2212
 394       q(1,1)=0d0
 395       q(2,1)=0d0
 396       q(3,1)=ebeam*beta
 397       q(4,1)=ebeam
 398
 399       istup(1)=-1
 400       idup(1)=2212
 401       mothup(1,1)=0
 402       mothup(2,1)=0
 403       icolup(1,1)=0
 404       icolup(2,1)=0
 405       do j=1,4
 406          pup(j,1)=q(j,1)
 407       enddo
 408       pup(5,1)=dsqrt(q(4,1)**2-q(3,1)**2-q(2,1)**2-q(1,1)**2)
 409       vtimup(1)=0
 410       spinup(1)=9.
 411
 412       q(1,2)=0d0
 413       q(2,2)=0d0
 414       q(3,2)=-ebeam*beta
 415       q(4,2)=ebeam
 416
 417       istup(2)=-1
 418       if(ppbar.eq.'true')then
 419          idup(2)=-2212
 420       else
 421          idup(2)=2212
 422       endif
 423       mothup(1,2)=0
 424       mothup(2,2)=0
 425       icolup(1,2)=0
 426       icolup(2,2)=0
 427       do j=1,4
 428          pup(j,2)=q(j,2)
 429       enddo
 430       pup(5,2)=dsqrt(q(4,2)**2-q(3,2)**2-q(2,2)**2-q(1,2)**2)
 431       vtimup(2)=0
 432       spinup(2)=9.
 433
 434 cccc  outgoing initial info
 435
 436       istup(3)=1
 437       if(ppbar.eq.'true')then
 438          idup(2)=-2212
 439       else
 440          idup(2)=2212
 441       endif
 442       idup(3)=2212
 443       mothup(1,3)=1
 444       mothup(2,3)=2
 445       icolup(1,3)=0
 446       icolup(2,3)=0
 447       vtimup(3)=0
 448       spinup(3)=9
 449
 450       istup(4)=1
 451       idup(4)=2212
 452       mothup(1,4)=1
 453       mothup(2,4)=2
 454       icolup(1,4)=0
 455       icolup(2,4)=0
 456       vtimup(4)=0
 457       spinup(4)=9
 458
 459 ccc   HEPEVT
 460       if(output.eq.'hepevt')then
 461
 462       nhep=nup
 463        write(6,*) "hepevt", nhep, nup
 464       do k=1,5
 465          phep(k,1)=pup(k,1)
 466          phep(k,2)=pup(k,2)
 467       enddo
 468
 469       do k=1,nup
 470          isthep(k)=istup(k)
 471          idhep(k)=idup(k)
 472          jmohep(1,k)=mothup(1,k)
 473          jmohep(2,k)=jmohep(1,k)
 474          jdahep(1,k)=0
 475          jdahep(2,k)=0
 476          vhep(1,k)=0d0
 477          vhep(2,k)=0d0
 478          vhep(3,k)=0d0
 479          vhep(4,k)=0d0
 480       enddo
 481
 482       if(pflag.eq.'rho')then
 483       jdahep(1,5)=7
 484       jdahep(2,5)=8
 485       jdahep(1,6)=9
 486       jdahep(2,6)=10
 487       endif
 488
 489       jmohep(2,5)=2
 490       jmohep(2,6)=2
 491
 492       endif
 493
 494       if(unw.eq.'true')then
 495          lmax=2
 496       else
 497          lmax=1
 498       endif
 499
 500 c      do ll=1,lmax
 501
 502       if(ll.eq.2)then
 503 c         ntotal=nev*10
 504          ntotal=1000000000
 505       endif
 506
 507       ip=ntotal+1
 508
 509 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 510 c                                                       c
 511 c     Start of event loop                               c
 512 c                                                       c
 513 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 514
 515       if(ll.eq.1)then
 516          write(6,*)'Generating weighted events...'
 517       else
 518         write(6,*)'Generating unweighted events...'
 519       endif
 520       return
 521       end
 522
 523       subroutine dimegenerate
 524       implicit none
 525
 526 ccccc hepevt output
 527       integer nmxhep,kk
 528       parameter (nmxhep=4000)
 529       integer nevhep,nhep,isthep,idhep,jmohep,jdahep
 530       double precision phep,vhep
 531       common /hepevt/ nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
 532      &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
 533 ccccc Les Houches  Event Common Block
 534       INTEGER MAXNUP
 535       PARAMETER (MAXNUP=500)
 536       INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP
 537       DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP
 538       COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
 539      &              IDUP(MAXNUP),ISTUP(MAXNUP),MOTHUP(2,MAXNUP),
 540      &              ICOLUP(2,MAXNUP),PUP(5,MAXNUP),VTIMUP(MAXNUP),
 541      &              SPINUP(MAXNUP)
 542 ccccc rambo variables
 543       integer npart,nrun
 544       double precision pin(4), am(100), pout(4,100), wt, ein
 545       common/ramboc/ pin, am, pout, wt, ein,
 546      &npart, nrun
 547
 548 ccccc
 549       double precision etaelmax,etaelmin,ptelmin,ptphmin,ecut,rmax,
 550      &rmin,mcut
 551       common/cuts/etaelmax,etaelmin,ptelmin,ptphmin,ecut,rmax,rmin,mcut
 552 ccccc
 553       double precision q(4,20)
 554       common /mom/ q
 555 ccccc
 556       double precision pt1(2), pt2(2), ptx(2)
 557       common/mompt/pt1,pt2,ptx
 558 ccccc
 559       double precision s, rts, mmes, yx
 560       double precision bjac, bp
 561       double precision mf127, mf1525
 562
 563       double precision m0, mmes0, mmes1, mmes2, mp, mwidth, pi, rt2,
 564      &ebeam, sh, th, uh, sum, sum1, weightm
 565
 566       integer id(20), ncut, ll, icut, nev, num, iin
 567
 568       common/vars/s,rts,mmes,yx
 569       common/hvars/sh,th,uh
 570       common/vars0/ mf127, mf1525, m0, mmes0, mmes1, mmes2, mp,
 571      &              mwidth, pi, rt2, ebeam, sum, sum1, weightm,
 572      &              bjac, bp
 573       common /ivars/ ncut, ll, icut, nev, num
 574 ccccc
 575       character prefix*50,fsp*10,order*10,pflag*10,fsi*10,formf*10
 576      &,ppbar*10,output*10,mregge*10,cuts*10,unw*10
 577       common/flags/ iin, pflag, fsi, ppbar, output, cuts, unw
 578       common/ff/formf
 579
 580 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 581 cccccc local variables ...
 582       double precision aa1, aa2, al1, al2, almax, almin, c, cc1, cc2,
 583      & msmax, msmin, cont, exn, mmesp, mwidth1, mwidth2, norm,
 584      & p1m, p1p, p2m, p2p, phi1, phi2, phix1, pt1sq, pt1x, pt1y,
 585      & pt2sq, pt2x, pt2y, ptxsq1, ptxsq2, ptxsqmax, ptxsqmin,
 586      & ran0, ran1, ran2, ran3, ran4, ran5, ran6, ranhist, ranx1,
 587      & ranx2, ranx3, ranx4, rm1, rm2, rmx1, rmx2,
 588      & root1sq, root2sq, snp, t1, t2, weight, weight0, weight1,
 589      & weight2, wthist, x1, x2, ymax1, ymax2, ymin1, ymin2, yx1,
 590      & yx2
 591       double precision pboo(4), pcm(4), plb(4)
 592       double precision svec(4)
 593       complex*16 zmat
 594       integer h, i, j, k, ntotal
 595       weight=0d0
 596       call r2455(ran0)
 597       call r2455(ran1)
 598       call r2455(ran2)
 599       call r2455(ran3)
 600       call r2455(ran4)
 601       call r2455(ran5)
 602       call r2455(ran6)
 603
 604 ccccccccccccccccccccccccccccccccccccccccccccccccccc
 605
 606 c     incoming protons
 607
 608       ID(1)=2212
 609       q(1,1)=0d0
 610       q(2,1)=0d0
 611       q(3,1)=ebeam
 612       q(4,1)=ebeam
 613
 614       ID(2)=2212
 615       q(1,2)=0d0
 616       q(2,2)=0d0
 617       q(3,2)=-ebeam
 618       q(4,2)=ebeam
 619
 620       phi1=2d0*pi*ran0
 621       phi2=2d0*pi*ran1
 622
 623       pt1sq=-dlog(ran2)/(bjac)
 624       pt2sq=-dlog(ran3)/(bjac)
 625
 626       weight=(dexp(bjac*pt1sq)*dexp(bjac*pt2sq))/bjac**2
 627
 628       pt1(1)=dsqrt(pt1sq)*dsin(phi1)
 629       pt1(2)=dsqrt(pt1sq)*dcos(phi1)
 630       pt2(1)=dsqrt(pt2sq)*dsin(phi2)
 631       pt2(2)=dsqrt(pt2sq)*dcos(phi2)
 632
 633       pt1x=pt1(1)
 634       pt1y=pt1(2)
 635       pt2x=pt2(1)
 636       pt2y=pt2(2)
 637
 638 cccccccccccccccccccccccccccccccc
 639
 640 ccc non-zero rho width
 641
 642 cccccccccccccccccccccccccccccccc
 643
 644        if(pflag.eq.'rho')then  ! new part here
 645
 646  677     call r2455(rm1)
 647          call r2455(rm2)
 648
 649          msmax=mmes0+4d0*mwidth
 650          msmin=2d0*0.13957018d0
 651
 652 cccccccccccccccc
 653
 654          almin=datan(-(mmes0**2-msmin**2)/mwidth/mmes0)
 655          almax=datan(-(mmes0**2-msmax**2)/mwidth/mmes0)
 656
 657          al1=almin+(almax-almin)*rm1
 658          al2=almin+(almax-almin)*rm2
 659
 660          mmes1=dsqrt(dtan(al1)*mmes0*mwidth+mmes0**2)
 661          mmes2=dsqrt(dtan(al2)*mmes0*mwidth+mmes0**2)
 662
 663          weight=weight*(almax-almin)**2
 664          weight=weight*mwidth**2*mmes0**2
 665          weight=weight*(1d0+dtan(al1)**2)*(1d0+dtan(al2)**2)
 666          weight=weight/4d0/mmes1/mmes2
 667
 668 ccccccccccccccccc
 669
 670          mwidth1=mwidth*((1d0-4d0*0.13957018d0**2/mmes1**2)/
 671      &        (1d0-4d0*0.13957018d0**2/mmes0**2))**1.5d0
 672          mwidth2=mwidth*((1d0-4d0*0.13957018d0**2/mmes2**2)/
 673      &        (1d0-4d0*0.13957018d0**2/mmes0**2))**1.5d0
 674
 675          if(mmes1.lt.2d0*0.13957018d0) goto 677
 676          if(mmes2.lt.2d0*0.13957018d0) goto 677
 677
 678          weight=weight*2d0*mmes0*mmes1*mwidth1/pi
 679          weight=weight*2d0*mmes0*mmes2*mwidth2/pi
 680          weight=weight*mmes1**2*mmes2**2/mmes0**4
 681          weight=weight/((mmes0**2-mmes1**2)**2+mmes1**4*
 682      &        mwidth1**2/mmes0**2)
 683          weight=weight/((mmes0**2-mmes2**2)**2+mmes2**4*
 684      &        mwidth2**2/mmes0**2)
 685
 686       endif
 687
 688
 689 cccccccccccccccccccccccccccccccccccccccccccccccccccc
 690
 691 c     pi rapidities
 692
 693          call r2455(ranx1)
 694          call r2455(ranx2)
 695          call r2455(ranx3)
 696          call r2455(ranx4)
 697
 698          phix1=2d0*pi*ranx1
 699
 700 c         r2max=1d0
 701 c         r2=r2max*ranx1
 702 c         ptxsq1=-dlog(r2)/bjac1  ! generate exp p_t^2 (can be more efficient)
 703
 704          ptxsqmin=0d0
 705          ptxsqmax=10d0         ! generate linear p_t^2
 706
 707          ptxsq1=ptxsqmin+(ptxsqmax-ptxsqmin)*ranx2
 708
 709          q(1,5)=dsqrt(ptxsq1)*dcos(phix1)
 710          q(2,5)=dsqrt(ptxsq1)*dsin(phix1)
 711          q(1,6)=-pt1(1)-pt2(1)-q(1,5)
 712          q(2,6)=-pt1(2)-pt2(2)-q(2,5)
 713          ptxsq2=q(1,6)**2+q(2,6)**2
 714
 715          rmx1=dsqrt(mmes1**2+ptxsq1)
 716          rmx2=dsqrt(mmes2**2+ptxsq2)
 717
 718          ymax1=dlog(rts/rmx1)
 719
 720          if(cuts.eq.'true')then
 721             if(rmax.lt.ymax1)then
 722                ymax1=rmax
 723             endif
 724          endif
 725
 726          ymin1=-ymax1
 727
 728          yx1=ymin1+(ymax1-ymin1)*ranx3
 729
 730          ymax2=(dlog((rts-rmx1*dexp(yx1))/rmx2))
 731          ymin2=-(dlog((rts-rmx1*dexp(-yx1))/rmx2))
 732
 733          if(cuts.eq.'true')then
 734             if(ymax2.gt.rmax)then
 735                ymax2=rmax
 736             endif
 737             if(ymin2.lt.-rmax)then
 738                ymin2=-rmax
 739             endif
 740          endif
 741
 742          if(ymax2.lt.ymin2)then
 743             weight=0d0
 744             goto 700
 745          endif
 746
 747          yx2=ymin2+(ymax2-ymin2)*ranx4
 748          x1=(rmx2*dexp(yx2)+rmx1*dexp(yx1))/rts
 749          x2=(rmx2*dexp(-yx2)+rmx1*dexp(-yx1))/rts
 750
 751          weight=weight*(ptxsqmax-ptxsqmin)
 752 c         weight=weight*dexp(bjac1*ptxsq1)*r2max/bjac1
 753          weight=weight*(ymax1-ymin1)*(ymax2-ymin2)
 754
 755
 756          q(3,5)=rmx1*(dexp(yx1)-dexp(-yx1))/2d0
 757          q(4,5)=rmx1*(dexp(yx1)+dexp(-yx1))/2d0
 758
 759          q(3,6)=rmx2*(dexp(yx2)-dexp(-yx2))/2d0
 760          q(4,6)=rmx2*(dexp(yx2)+dexp(-yx2))/2d0
 761
 762 ccccccccccccccccccccccccccccccccccccccccccccccccccccc
 763
 764 c     impose massive on-shell condition by solving
 765 c                   p1+ + cc1/p2- = aa1
 766 c                   p2- + cc2/p1+ = aa2
 767
 768       aa1=bp*(1d0-x1)
 769       aa2=bp*(1d0-x2)
 770       cc1=0.5d0*(pt2sq+mp**2)
 771       cc2=0.5d0*(pt1sq+mp**2)
 772
 773       root1sq=(cc1-cc2-aa1*aa2)**2-4d0*cc2*aa1*aa2
 774       root2sq=(cc2-cc1-aa1*aa2)**2-4d0*cc1*aa1*aa2
 775       if(root1sq.le.0d0.or.root2sq.le.0d0)then
 776          weight=0d0
 777          goto 700
 778       endif
 779       p1p=(cc2-cc1+aa1*aa2+dsqrt(root1sq))/(2d0*aa2)
 780       p2m=(cc1-cc2+aa1*aa2+dsqrt(root2sq))/(2d0*aa1)
 781       p1m=(pt1sq+mp**2)/(2d0*p1p)
 782       p2p=(pt2sq+mp**2)/(2d0*p2m)
 783
 784       if(p1p.lt.0d0.or.p1m.lt.0d0.or.p2p.lt.0d0.or.p2m.lt.0d0)then
 785          weight=0d0
 786          goto 700
 787       endif
 788
 789       t1=-rts*p1m*rt2
 790       t2=-rts*p2p*rt2
 791
 792       q(1,3)=pt1(1)
 793       q(2,3)=pt1(2)
 794       q(3,3)=(p1p-p1m)/rt2
 795       q(4,3)=(p1p+p1m)/rt2
 796
 797       q(1,4)=pt2(1)
 798       q(2,4)=pt2(2)
 799       q(3,4)=(p2p-p2m)/rt2
 800       q(4,4)=(p2p+p2m)/rt2
 801
 802 ccccccccccccccccccccccccccccccccccccccccccccccccc
 803
 804       do k=1,4
 805          svec(k)=q(k,5)+q(k,6)
 806       enddo
 807
 808          call dot(svec,svec,sh)
 809
 810 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 811
 812 c     rho0 decays
 813
 814 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 815
 816          if(pflag.eq.'rho')then
 817
 818             kk=6
 819
 820             do k=5,6
 821
 822                if(k.eq.5)then
 823                   mmesp=mmes1
 824                else
 825                   mmesp=mmes2
 826                endif
 827                call rambo(npart,mmesp,am,pout,wt) ! call RAMBO
 828
 829                do j=1,4
 830                   pboo(j)=q(j,k)
 831                enddo
 832
 833                do h=1,2
 834                   do j=1,4
 835                      pcm(j)=pout(j,h)
 836                   enddo
 837                   call boost(mmesp,pboo,pcm,plb)
 838                   kk=kk+1
 839                   do j=1,4
 840                      q(j,kk)=plb(j)
 841                   enddo
 842                enddo
 843
 844             enddo
 845
 846          endif
 847
 848
 849 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 850
 851 c      place cuts
 852
 853          if(cuts.eq.'true')then
 854
 855             call cut(icut)
 856
 857             if(icut.eq.0)then
 858                weight=0d0
 859                weight0=0d0
 860                weight1=0d0
 861                weight2=0d0
 862                ncut=ncut+1
 863                goto 700
 864             endif
 865
 866          endif
 867
 868 c         print*,i
 869
 870 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 871
 872 c     Write 4-momenta to event record array
 873
 874 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 875
 876          if(output.eq.'lhe')then
 877
 878          do k=3,4
 879             do j=1,4
 880                pup(j,k)=q(j,k)
 881             enddo
 882             pup(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2)
 883          enddo
 884
 885          if(pflag.eq.'rho')then
 886
 887             do k=5,10
 888                do j=1,4
 889                   pup(j,k)=q(j,k)
 890                enddo
 891                pup(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2)
 892             enddo
 893
 894          else
 895
 896             do k=5,6
 897                do j=1,4
 898                   pup(j,k)=q(j,k)
 899                enddo
 900                pup(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2)
 901             enddo
 902
 903          endif
 904
 905          else
 906
 907          do k=3,4
 908             do j=1,4
 909                phep(j,k)=q(j,k)
 910             enddo
 911             phep(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2)
 912          enddo
 913
 914          if(pflag.eq.'rho')then
 915
 916             do k=5,10
 917                do j=1,4
 918                   phep(j,k)=q(j,k)
 919                enddo
 920                phep(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2)
 921             enddo
 922
 923          else
 924
 925             do k=5,6
 926                do j=1,4
 927                   phep(j,k)=q(j,k)
 928                enddo
 929                phep(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2)
 930             enddo
 931
 932          endif
 933
 934          endif
 935
 936 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 937
 938 c     Weight evaluation
 939
 940 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 941
 942 c     pp-->pp(M_1 M_2) matrix element
 943
 944 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 945
 946
 947       call schimc(pt1x,pt1y,pt2x,pt2y,zmat)
 948
 949       weight=weight*cdabs(zmat)**2
 950       weight=weight/norm**2
 951
 952 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 953
 954       weight=weight/(s**2*(16d0)**3*pi**5)*389d3
 955
 956 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 957
 958 ccc   Probability of no additional particles produced in production subprocess
 959
 960 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 961
 962       if(fsi.eq.'true')then
 963
 964          if(pflag.eq.'pipm'.or.pflag.eq.'pi0')then
 965
 966             if(dsqrt(sh).gt.mf127)then
 967                m0=mf127
 968                c=0.7d0
 969                cont=1d0-(m0-2d0*mmes)**2/m0**2
 970                exn=c*dlog(((dsqrt(sh)-2d0*mmes)**2)/m0**2+cont)
 971                snp=dexp(-exn)
 972             else
 973                snp=1d0
 974             endif
 975
 976          elseif(pflag.eq.'kpkm'.or.pflag.eq.'ks')then
 977
 978             if(dsqrt(sh).gt.mf1525)then
 979                m0=mf1525
 980                c=0.7d0
 981                cont=1d0-(m0-2d0*mmes)**2/m0**2
 982                exn=c*dlog(((dsqrt(sh)-2d0*mmes)**2)/m0**2+cont)
 983                snp=dexp(-exn)
 984             else
 985                snp=1d0
 986             endif
 987
 988          else
 989
 990             snp=1d0
 991
 992          endif
 993
 994       else
 995
 996          snp=1d0
 997
 998       endif
 999
1000       if(snp.gt.1d0)then
1001          snp=1d0
1002       endif
1003
1004       weight=weight*snp
1005
1006 ccccccccccccccc
1007
1008       if(pflag.eq.'pi0'.or.pflag.eq.'ks'.or.pflag.eq.'rho')then
1009          weight=weight/2d0    ! symmetry factor
1010       endif
1011
1012 cccccccccccccccccccccccccccccccccccccccccccccccccccccc
1013
1014  666  weight=weight
1015
1016       if(weightm.lt.weight)then
1017          weightm=weight
1018       endif
1019
1020
1021       if(ll.eq.1)then
1022
1023          xwgtup=weight
1024
1025          wthist=weight/dfloat(ntotal)
1026
1027          call binit(wthist)
1028
1029       else
1030
1031          call r2455(ranhist)
1032
1033          if(ranhist.lt.weight/weightm)then
1034
1035             xwgtup=1d0
1036
1037             num=num+1
1038
1039
1040 c            call binit(sumt/nev)
1041
1042             write(35,302)num,nup
1043
1044             if(output.eq.'lhe')then
1045
1046             do j=1,nup
1047                write(35,300)j,istup(j),idup(j),mothup(1,j),mothup(2,j),
1048      &              icolup(1,j),icolup(2,j),pup(1,j),pup(2,j),
1049      &              pup(3,j),pup(4,j),pup(5,j),vtimup(j),spinup(j)
1050             enddo
1051
1052             else
1053
1054                do j=1,nup
1055                   write(35,301)j,idhep(j),isthep(j),jmohep(1,j),
1056      &                 jmohep(2,j),jdahep(1,j),jdahep(2,j),
1057      &                 phep(1,j),phep(2,j),phep(3,j),phep(4,j),phep(5,j)
1058      &                 ,vhep(1,j),vhep(2,j),vhep(3,j),vhep(4,j)
1059                enddo
1060
1061             endif
1062
1063
1064          endif
1065
1066       endif
1067
1068       if(ll.eq.2)then
1069
1070       if(num.gt.nev-1)then  ! exit loop once required no. of unweighted events generated
1071
1072          ntotal=i
1073 !         goto 888
1074          call terminate()
1075
1076       endif
1077       endif
1078
1079  700  sum=sum+weight
1080       sum1=sum1+weight**2
1081
1082 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1083 c                                                                    c
1084 c     End of event loop                                              c
1085 c                                                                    c
1086 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1087       call terminate()
1088  300  format(i4,1x,i4,1x,i8,1x,i4,1x,i4,1x,i4,1x,i4,1x,E13.6,1x,
1089      &E13.6,1x,E13.6,1x,E13.6,1x,E13.6,1x,E13.6,1x,E13.6)
1090  301  format(i4,1x,i8,1x,i4,1x,i4,1x,i4,1x,i4,1x,i4,1x,E13.6,1x,E13.6
1091      &,1x,E13.6,1x,E13.6,1x,E13.6,1x,E13.6,1x,E13.6,1x,E13.6,1x
1092      &,E13.6)
1093  302  format(1x,i8,1x,i5,1x,E13.6)
1094
1095       return
1096       end
1097
1098       subroutine terminate
1099       character prefix*50,fsp*10,order*10,pflag*10,fsi*10,formf*10
1100      &,ppbar*10,output*10,mregge*10,cuts*10,unw*10
1101
1102  888  continue
1103       sum=sum/dfloat(ntotal)
1104       sum1=sum1/dfloat(ntotal)
1105       var=dsqrt((sum1-sum**2)/dfloat(ntotal))
1106       EFF=1d0*(ntotal-ncut)/ntotal
1107
1108       if(ll.eq.1)then
1109          print*,'...done!'
1110          write(6,*)
1111          sumt=sum
1112       else
1113          print*,'...done!'
1114          write(6,*)
1115       endif
1116
1117
1118
1119 ccccc     create histograms
1120
1121       if(ll.eq.1)then
1122
1123       do i=1,nhist
1124            call histo2(i,0)
1125       enddo
1126
1127       endif
1128
1129       if(ll.eq.1)then
1130       if(pflag.eq.'pipm')then
1131          write(6,*)
1132          write(6,*)'pi+pi- production'
1133       elseif(pflag.eq.'pi0')then
1134          write(6,*)
1135          write(6,*)'pi0pi0 production'
1136       elseif(pflag.eq.'kpkm')then
1137          write(6,*)
1138          write(6,*)'K+K- production'
1139       elseif(pflag.eq.'ks')then
1140          write(6,*)
1141          write(6,*)'K0K0 production'
1142       elseif(pflag.eq.'rho')then
1143          write(6,*)
1144          write(6,*)'rho0rho0 production'
1145       endif
1146       endif
1147
1148       if(ll.eq.1)then
1149          if(pflag.eq.'rho')then
1150          write(6,221)rts,ntotal,sum,var
1151       else
1152          write(6,222)rts,mmes,ntotal,sum,var
1153       endif
1154       else
1155          if(output.eq.'lhe')then
1156             print*,'LHE output'
1157          else
1158             print*,'HEPEVT output'
1159          endif
1160          write(6,223)nev,sum,var
1161       endif
1162
1163 !      enddo
1164
1165       close(35)
1166
1167  221  format(/,
1168      .       3x,'collider energy  (GeV)    : ',f10.3,/,
1169      .       3x,'number of events          : ',i9,/,
1170      .       3x,'sigma (nb)                : ',f15.5,'  +-  ',f15.5,
1171      .   /)
1172
1173  222  format(/,
1174      .       3x,'collider energy  (GeV)    : ',f10.3,/,
1175      .       3x,'meson mass (GeV)          : ',f10.5,/,
1176      .       3x,'number of events          : ',i9,/,
1177      .       3x,'sigma (nb)                : ',f15.5,'  +-  ',f15.5,
1178      .   /)
1179
1180  223  format(3x,'number of events          : ',i9,/,
1181      .       3x,'sigma (nb)                : ',f15.5,'  +-  ',f15.5,
1182      .   /)
1183
1184       call cpu_time(t2)
1185       print*,'time elapsed = ', t2, ' s'
1186
1187       return
1188       end
1189
1190 ccccc Pomeron -- (off-shell) meson form factor
1191
1192       function fpi(x)
1193       implicit double precision(a-y)
1194       character formf*10
1195       common/vars/s,rts,mmes,yx
1196       common/ff/formf
1197       common/ffpars/bexp,ao,bo,aoo,boo
1198
1199       if(formf.eq.'exp')then
1200          fpi=exp((x-mmes**2)*bexp)
1201       elseif(formf.eq.'orexp')then
1202          fpi=exp(-(dsqrt(-x+mmes**2+ao**2)-ao)*bo)
1203       elseif(formf.eq.'power')then
1204          fpi=1d0/(1d0-(x-mmes**2)/aoo)
1205       endif
1206
1207       return
1208       end
1209
1210 ccccccc  Pom Pom --> meson pair amplitude
1211
1212       subroutine wev(matt,v1,v2)
1213       implicit double precision(a-y)
1214       implicit complex*16(z)
1215       double precision q(4,20),svec(4),tvec(4),uvec(4),v1(4),v2(4)
1216       integer k
1217       character*10 mregge,pflag
1218       complex*16 matt
1219       common/mom/q
1220       common/vars/s,rts,mmes,yx
1221       common/wvars/sig0,bb,t1,t2
1222       common/alphas/alphap,alpha0,alphapr,alpha0r,alphapm,alpha0m
1223       common/sec/cpom,cf,crho,aff,ar
1224       common/hvars/sh,th,uh
1225       common/regge/mregge
1226       common/flags/ iin, pflag, fsi, ppbar, output, cuts, unw
1227
1228
1229       zi=(0d0,1d0)
1230       pi=dacos(-1d0)
1231
1232       do k=1,4
1233          q(k,10)=v1(k)
1234          q(k,11)=v2(k)
1235       enddo
1236
1237       do k=1,4
1238          tvec(k)=q(k,10)-q(k,5)
1239          uvec(k)=q(k,10)-q(k,6)
1240       enddo
1241
1242          call dot(tvec,tvec,th)
1243          call dot(uvec,uvec,uh)
1244
1245          sht1=2d0*(q(4,5)*q(4,3)-q(3,5)*q(3,3)-q(2,5)*q(2,3)-
1246      &q(1,5)*q(1,3))+mmes**2
1247
1248          sht2=2d0*(q(4,6)*q(4,4)-q(3,6)*q(3,4)-q(2,6)*q(2,4)-
1249      &q(1,6)*q(1,4))+mmes**2
1250
1251          shu1=2d0*(q(4,6)*q(4,3)-q(3,6)*q(3,3)-q(2,6)*q(2,3)-
1252      &q(1,6)*q(1,3))+mmes**2
1253
1254          shu2=2d0*(q(4,5)*q(4,4)-q(3,5)*q(3,4)-q(2,5)*q(2,4)-
1255      &q(1,5)*q(1,4))+mmes**2
1256
1257        bt1=(2d0*alphap*dlog(sht1))
1258        bt2=(2d0*alphap*dlog(sht2))
1259        bu1=(2d0*alphap*dlog(shu1))
1260        bu2=(2d0*alphap*dlog(shu2))
1261
1262        bt1r=(2d0*alphapr*dlog(sht1))
1263        bt2r=(2d0*alphapr*dlog(sht2))
1264        bu1r=(2d0*alphapr*dlog(shu1))
1265        bu2r=(2d0*alphapr*dlog(shu2))
1266
1267        sig0t1=sig0*sht1**0.0808d0/0.389d0
1268        sig0t2=sig0*sht2**0.0808d0/0.389d0
1269        sig0u1=sig0*shu1**0.0808d0/0.389d0
1270        sig0u2=sig0*shu2**0.0808d0/0.389d0
1271
1272        t11=v1(1)**2+v1(2)**2
1273        t22=v2(1)**2+v2(2)**2
1274
1275           zmt=(zi*dexp(-bt1*t11/2d0)*cpom*sht1**alpha0
1276      &         +((aff+zi)*cf*sht1**alpha0r+(ar-zi)*crho*sht1**alpha0r)
1277      &         *dexp(-bt1r*t11/2d0))
1278           zmt=zmt*(zi*dexp(-bt2*t22/2d0)*cpom*sht2**alpha0
1279      &         +((aff+zi)*cf*sht2**alpha0r-(ar-zi)*crho*sht2**alpha0r)
1280      &         *dexp(-bt2r*t11/2d0))
1281           zmt=zmt*fpi(th)**2/(mmes**2-th)
1282
1283           zmu=(zi*dexp(-bu1*t11/2d0)*cpom*shu1**alpha0
1284      &         +((aff+zi)*cf*shu1**alpha0r-(ar-zi)*crho*shu1**alpha0r)
1285      &         *dexp(-bu1r*t11/2d0))
1286           zmu=zmu*(zi*dexp(-bu2*t22/2d0)*cpom*shu2**alpha0
1287      &         +((aff+zi)*cf*shu2**alpha0r+(ar-zi)*crho*shu2**alpha0r)
1288      &         *dexp(-bu2r*t22/2d0))
1289           zmu=zmu*fpi(uh)**2/(mmes**2-uh)
1290
1291           matt=zmu+zmt
1292
1293        return
1294        end
1295
1296 c     binning subroutine
1297
1298       subroutine binit(wt)
1299       implicit double precision(a-y)
1300       double precision q(4,20),pt1(2),pt2(2),ptx(2)
1301       common/mom/q
1302       common/mompt/pt1,pt2,ptx
1303       common/vars/s,rts,mmes,yx
1304       common/vars1/ptgam,etagam,ptel2,ptel1,etael1,etael2
1305       common/vars2/ptpi1,etapi1,ptpi2,etapi2
1306       common/ang/cost,costa,costb
1307       common/hvars/sh,th,uh
1308
1309       ypisum=0.5d0*dlog((q(4,5)+q(4,6)+q(3,5)+q(3,6))
1310      &     /(q(4,5)+q(4,6)-q(3,5)-q(3,6)))
1311
1312       ypi1=0.5d0*dlog((q(4,5)+q(3,5))
1313      &     /(q(4,5)-q(3,5)))
1314       ypi2=0.5d0*dlog((q(4,6)+q(3,6))
1315      &     /(q(4,6)-q(3,6)))
1316
1317       call histo1(1,30,0d0,4.5d0,dsqrt(sh),wt)
1318
1319       return
1320       end
1321
1322 ccccc Cutting subroutine
1323 ccccc
1324 ccccc User can define their own cuts on any particle momenta, stored in array q(i,j)
1325 ccccc
1326
1327       subroutine cut(icut)
1328       implicit double precision(a-y)
1329       double precision q(4,20)
1330       integer icut
1331       common/mom/q
1332       common/vars/s,rts,mmes,yx
1333       common/vars1/ptgam,etagam,ptel2,ptel1,etael1,etael2
1334       common/vars2/ptpi1,etapi1,ptpi2,etapi2
1335       common/cuts/etaelmax,etaelmin,ptelmin,ptphmin,ecut,rmax,rmin,mcut
1336       common/hvars/sh,th,uh
1337       integer iii
1338       common/ivar/iii
1339
1340       icut=0
1341
1342 c -- meson 1 variables
1343       ptpi1=dsqrt(q(1,5)**2+q(2,5)**2)
1344       pmodpi1=dsqrt(q(1,5)**2+q(2,5)**2+q(3,5)**2)
1345       etapi1=.5d0*dlog((pmodpi1+q(3,5))
1346      .              /(pmodpi1-q(3,5)))
1347
1348       ypi1=.5d0*dlog((q(4,5)+q(3,5))
1349      .              /(q(4,5)-q(3,5)))
1350
1351         etpi1=q(4,5)*dsqrt(pmodpi1**2-q(3,5)**2)/pmodpi1
1352
1353 c -- meson 2 variables
1354       ptpi2=dsqrt(q(1,6)**2+q(2,6)**2)
1355       pmodpi2=dsqrt(q(1,6)**2+q(2,6)**2+q(3,6)**2)
1356       etapi2=.5d0*dlog((pmodpi2+q(3,6))
1357      .              /(pmodpi2-q(3,6)))
1358
1359       ypi2=.5d0*dlog((q(4,6)+q(3,6))
1360      .              /(q(4,6)-q(3,6)))
1361
1362         etpi2=q(4,6)*dsqrt(pmodpi2**2-q(3,6)**2)/pmodpi2
1363
1364 c        print*,ptpi1,ptpi2
1365
1366 ccccc    example cuts....
1367
1368 c        if(dsqrt(sh).lt.0.8d0)return
1369
1370         xf1=2d0*q(3,3)/rts
1371         xf2=2d0*q(3,4)/rts
1372
1373 c        if(dabs(xf1).lt.0.9d0)return
1374 c        if(dabs(xf2).lt.0.9d0)return
1375
1376 ccccccccc  Basic cuts described at start of code
1377
1378         if(ptpi1.lt.ecut)return
1379         if(ptpi2.lt.ecut)return
1380         if(ypi1.gt.rmax)return
1381         if(ypi2.gt.rmax)return
1382         if(ypi1.lt.rmin)return
1383         if(ypi2.lt.rmin)return
1384
1385         icut=1
1386
1387       return
1388       end
1389
1390 c     prints histograms
1391
1392       subroutine histo1(ih,ib,x0,x1,x,w)
1393       implicit real*8(a-h,o-z)
1394       character*1 regel(30),blank,star
1395       dimension h(20,100),hx(20),io(20),iu(20),ii(20)
1396       dimension y0(20),y1(20),ic(20)
1397       data regel / 30*' ' /,blank /' ' /,star /'*'/
1398       save
1399       open(10,file="output.dat")
1400       y0(ih)=x0
1401       y1(ih)=x1
1402       ic(ih)=ib
1403       if(x.lt.x0) goto 11
1404       if(x.gt.x1) goto 12
1405       ix=idint((x-x0)/(x1-x0)*dble(ib))+1
1406       h(ih,ix)=h(ih,ix)+w
1407       if(h(ih,ix).gt.hx(ih)) hx(ih)=h(ih,ix)
1408       ii(ih)=ii(ih)+1
1409       return
1410    11 iu(ih)=iu(ih)+1
1411       return
1412    12 io(ih)=io(ih)+1
1413       return
1414       entry histo2(ih,il)
1415       ib1=ic(ih)
1416       x01=y0(ih)
1417       x11=y1(ih)
1418       bsize=(x11-x01)/dble(ib1)
1419       hx(ih)=hx(ih)*(1.d0+1.d-06)
1420       if(il.eq.0) write(6,21) ih,ii(ih),iu(ih),io(ih)
1421       if(il.eq.1) write(6,22) ih,ii(ih),iu(ih),io(ih)
1422    21 format(' no.',i3,' lin : inside,under,over ',3i6)
1423    22 format(' no.',i3,' log : inside,under,over ',3i6)
1424       if(ii(ih).eq.0) goto 28
1425       write(6,23)
1426    23 format(35(1h ),3(10h----+----i))
1427       do 27 iv=1,ib1
1428       z=(dble(iv)-0.5d0)/dble(ib1)*(x11-x01)+x01
1429       if(il.eq.1) goto 24
1430       iz=idint(h(ih,iv)/hx(ih)*30.)+1
1431       goto 25
1432    24 iz=-1
1433       if(h(ih,iv).gt.0.d0)
1434      .iz=idint(dlog(h(ih,iv))/dlog(hx(ih))*30.)+1
1435    25 if(iz.gt.0.and.iz.le.30) regel(iz)=star
1436       write(6,26) z,h(ih,iv)/bsize,(regel(i),i=1,30)
1437 c      write(10,*)z-0.125d0/2d0,h(ih,iv)/bsize        ! Print histogram to file
1438       write(10,*)z,h(ih,iv)/bsize        ! Print histogram to file
1439 c      write(10,*)z,h(ih,iv)/bsize        ! Print histogram to file
1440    26 format(1h ,2g15.6,4h   i,30a1,1hi)
1441    36 format(1h ,2g15.6)
1442       if(iz.gt.0.and.iz.le.30) regel(iz)=blank
1443    27 continue
1444       write(6,23)
1445       return
1446    28 write(6,29)
1447    29 format('  no entries inside histogram')
1448       return
1449       entry histo3(ih)
1450       do 31 i=1,100
1451    31 h(ih,i)=0.
1452       hx(ih)=0.
1453       io(ih)=0
1454       iu(ih)=0
1455       ii(ih)=0
1456       close(10)
1457       return
1458       end
1459
1460
1461 cccc  Initializes soft model parameters
1462
1463       subroutine initpars(in)
1464       implicit double precision(a-y)
1465       integer in,i1,i2
1466       integer nch
1467       double precision cc0(5),bm(5),bb0(5),pp0(5),bex(5),gaa(5)
1468       common/pars/cc0,bm,bb0,pp0,bex,asp,sigo,gaa,ep,norm
1469       common/ipars/nch
1470       common/vars/s,rts,mmes,yx
1471
1472       pi=dacos(-1d0)
1473
1474       if(in.eq.1)then
1475
1476          ep=0.13d0
1477          asp=0.08d0
1478          ep1=0d0
1479          sigo=23d0/0.39d0
1480          gaa(3)=0.4d0
1481          gd2=0.3d0
1482          nch=2
1483          ntf=0
1484
1485          cc0(1)=0.45d0
1486          bm(1)=3d0
1487          bb0(1)=0.1d0
1488          pp0(1)=0.92d0
1489          bex(1)=8.5d0
1490
1491          cc0(2)=0.45d0
1492          bm(2)=1.5d0
1493          bb0(2)=0.5d0
1494          pp0(2)=0.1d0
1495          bex(2)=4.5d0
1496
1497          cc0(3)=1d0
1498          bm(3)=0d0
1499          bb0(3)=0.8d0
1500          pp0(3)=0.5d0
1501          bex(3)=0.5d0
1502
1503       elseif(in.eq.2)then
1504
1505          ep=0.115d0
1506          asp=0.11d0
1507          ep1=0d0
1508          sigo=33d0/0.39d0
1509          gaa(3)=0.6d0
1510          gd2=0.16d0
1511          nch=2d0
1512          ntf=0d0
1513
1514          cc0(1)=0.63d0
1515          bm(1)=3d0
1516          bb0(1)=0.1d0
1517          pp0(1)=0.5d0
1518          bex(1)=8d0
1519
1520          cc0(2)=0.47d0
1521          bm(2)=1.5d0
1522          bb0(2)=0.5d0
1523          pp0(2)=0.1d0
1524          bex(2)=6d0
1525
1526          cc0(3)=1d0
1527          bm(3)=0d0
1528          bb0(3)=0.8d0
1529          pp0(3)=0.5d0
1530          bex(3)=0.5d0
1531
1532       elseif(in.eq.3)then
1533
1534          ep=0.093d0
1535          asp=0.075d0
1536          ep1=0d0
1537          sigo=60d0/0.39d0
1538          gaa3=4.8d0
1539          gd2=1.03d0
1540          nch=2
1541          nga=1
1542          cc0(1)=0.55d0
1543          bm(1)=3d0
1544          bb0(1)=0.27d0
1545          pp0(1)=0.48d0
1546          bex(1)=5.3d0
1547
1548          cc0(2)=0.48d0
1549          bm(2)=1.5d0
1550          bb0(2)=0.1d0
1551          pp0(2)=1d0
1552          bex(2)=3.8d0
1553
1554          cc0(3)=0.24d0
1555
1556       elseif(in.eq.4)then
1557
1558          ep=0.11d0              !!! Capital delta
1559          asp=0.06d0             !!! alpha'
1560          ep1=0d0                !!! Zero in all models, matters (?)
1561          sigo=50d0/0.39d0       !!! sigma_0 (GeV^-2)
1562          gaa3=6d0               !!! k2/k(1.8 TeV)
1563          gd2=1.3d0              !!! k1/k(1.8 TeV)
1564          nch=2
1565          nga=1                  !!! A flag (?)
1566          cc0(1)=0.6d0           !!! d1
1567          bm(1)=3d0              !!! Doesn't matter
1568          bb0(1)=0.45d0          !!! c1-0.08 (added back later)
1569          pp0(1)=0.5d0           !!! 2*|a_1|^2
1570          bex(1)=7.2d0           !!! b1
1571
1572          cc0(2)=0.48d0          !!! d2
1573          bm(2)=1.5d0            !!! Doesn't matter
1574          bb0(2)=0.16d0          !!! c2-0.08 (added back later)
1575          pp0(2)=1d0             !!! |a_2|^2 is set later
1576          bex(2)=4.2d0           !!! b2
1577
1578          cc0(3)=0.12d0          !!! Beta : k^2_min ~ s^Beta
1579
1580       endif
1581
1582       if(nch.eq.3) pp0(3)=3d0-pp0(2)-pp0(1)
1583       if(nch.eq.2) pp0(2)=2d0-pp0(1) !!! Set |a_2|^2
1584       if(nch.eq.1) pp0(1)=1d0
1585
1586       if(in.eq.3.or.in.eq.4)then
1587
1588       gamm=(1800d0/rts)**cc0(3)
1589       ga1=1d0/(1d0+gamm*gd2)
1590       ga2=1d0/(1d0+gamm*gaa3)
1591       gaa(1)=2d0*ga1/(ga1+ga2)   !!! gamma_1 (+) Eq (15)
1592       gaa(2)=2d0*ga2/(ga1+ga2)   !!! gamma_2 (-) Eq (15)
1593
1594       elseif(in.eq.1.or.in.eq.2)then
1595
1596          if(nch.eq.2)then
1597             gaa(1)=1d0+dsqrt(gd2)
1598             gaa(2)=1d0-dsqrt(gd2)
1599             gaa(3)=0.
1600          elseif(nch.eq.1)then
1601             gaa(1)=1d0
1602             gaa(2)=0d0
1603             gaa(3)=0d0
1604          endif
1605
1606       endif
1607
1608       sum=0d0
1609
1610       do i1=1,nch
1611          do i2=1,nch
1612             sum=sum+gaa(i1)*gaa(i2)*pp0(i1)*pp0(i2)/dble(nch)**2
1613          enddo
1614       enddo
1615
1616       norm=sum
1617
1618       return
1619       end
1620
1621 ccccc Calculates proton opacity
1622
1623       subroutine calcop
1624       implicit double precision(a-y)
1625       integer nb,i1,i2,nch,ib
1626       common/ipars/nch
1627       double precision op(5,5,10000,2),oph(5,5,10000,2)
1628       common/opac/op,oph
1629
1630       nb=900
1631       hb=100d0/dble(nb)
1632
1633       print*,'Calculating opacity'
1634
1635
1636       do ib=1,nb+1
1637          bt=dble(ib-1)*hb
1638
1639          do i1=1,nch
1640             do i2=1,nch
1641
1642                call opacity(i1,i2,bt,fr,fr1)
1643
1644
1645                op(i1,i2,ib,1)=bt
1646                op(i1,i2,ib,2)=fr
1647                oph(i1,i2,ib,1)=bt
1648                oph(i1,i2,ib,2)=fr1
1649
1650                write(40,*)bt,fr,fr1
1651
1652             enddo
1653          enddo
1654       enddo
1655
1656       return
1657       end
1658
1659 cccc  Calculates screening amplitude (in k_t space)
1660
1661       subroutine calcscreen
1662       implicit double precision(a-y)
1663       integer ns,i1,i2,nch,ib
1664       common/ipars/nch
1665       double precision sca(5,5,0:40000,2),sca1(5,5,0:40000,2)
1666       common/spac/sca,sca1
1667
1668       ns=900
1669       ksqma=8.2d0
1670       inck=ksqma/dble(ns)
1671       ksqmin=0.001d0
1672       lginck=(dlog(ksqma/ksqmin))/dble(ns)
1673
1674       print*,'Calculating screening amplitude'
1675
1676       do ib=0,ns+1
1677
1678          ksq=dble(ib-1)*inck
1679          lgksq=ksq
1680
1681          if(ib.eq.0)then
1682             ksq=0d0
1683             lgksq=0d0
1684          else
1685             lgksq=dble(ib-1)*lginck+dlog(ksqmin)
1686             ksq=dexp(lgksq)
1687          endif
1688
1689          do i1=1,nch
1690             do i2=1,nch
1691
1692                call screening(i1,i2,ksq,sc,sc1)
1693
1694                sca(i1,i2,ib,1)=lgksq
1695                sca(i1,i2,ib,2)=sc
1696                sca1(i1,i2,ib,1)=lgksq
1697                sca1(i1,i2,ib,2)=sc1
1698
1699             enddo
1700          enddo
1701
1702       enddo
1703
1704       return
1705       end
1706
1707 ccccc Interpolation....
1708
1709       subroutine opacityint(i,j,bt,out,out1)
1710       implicit double precision(a-y)
1711       double precision op(5,5,10000,2),oph(5,5,10000,2)
1712       integer i,j,it
1713       common/opac/op,oph
1714
1715       incbt=op(1,1,2,1)-op(1,1,1,1)
1716       it=nint(bt/incbt)
1717       if(dble(it).gt.(bt/incbt))then
1718          it=it-1
1719       endif
1720
1721       m=(op(i,j,it+2,2)-op(i,j,it+1,2))
1722      &/(op(i,j,it+2,1)-op(i,j,it+1,1))
1723       del=bt-op(1,1,it+1,1)
1724       mh=(oph(i,j,it+2,2)-oph(i,j,it+1,2))
1725      &/(oph(i,j,it+2,1)-oph(i,j,it+1,1))
1726       delh=bt-oph(1,1,it+1,1)
1727
1728       out=m*del+op(i,j,it+1,2)
1729       out1=mh*delh+oph(i,j,it+1,2)
1730
1731       return
1732       end
1733
1734       subroutine screeningint(i,j,ktsq,out,out1)
1735       implicit double precision(a-y)
1736       double precision sca(5,5,0:40000,2),sca1(5,5,0:40000,2)
1737       integer i,j,it,ns
1738       common/spac/sca,sca1
1739
1740       if(ktsq.lt.0.001d0)then
1741
1742          it=0
1743
1744          m=(sca(i,j,it+1,2)-sca(i,j,it,2))
1745      &        /(dexp(sca(i,j,it+1,1))-sca(i,j,it,1))
1746          del=ktsq-sca(1,1,it,1)
1747          m1=sca1(i,j,it+1,2)-sca1(i,j,it,2)
1748      &        /(dexp(sca1(i,j,it+1,1))-sca1(i,j,it,1))
1749          del1=ktsq-sca1(1,1,it,1)
1750
1751          out=m*del+sca(i,j,it,2)
1752          out1=m1*del1+sca1(i,j,it,2)
1753
1754       elseif(ktsq.lt.8d0)then
1755
1756          ktmin=sca(1,1,1,1)
1757          inckt=sca(1,1,2,1)-sca(1,1,1,1)
1758          it=nint((dlog(ktsq)-ktmin)/inckt)
1759          if(dble(it).gt.((dlog(ktsq)-ktmin)/inckt))then
1760             it=it-1
1761          endif
1762
1763          m=(sca(i,j,it+2,2)-sca(i,j,it+1,2))
1764      &        /(sca(i,j,it+2,1)-sca(i,j,it+1,1))
1765          del=dlog(ktsq)-sca(1,1,it+1,1)
1766          m1=sca1(i,j,it+2,2)-sca1(i,j,it+1,2)
1767      &        /(sca1(i,j,it+2,1)-sca1(i,j,it+1,1))
1768          del1=dlog(ktsq)-sca1(1,1,it+1,1)
1769
1770          out=m*del+sca(i,j,it+1,2)
1771          out1=m1*del1+sca1(i,j,it+1,2)
1772
1773       else
1774
1775          out=0d0
1776          out1=0d0
1777
1778       endif
1779
1780       return
1781       end
1782
1783 cccc  Integrates round Pomeron loop (to calculate screened amplitude)
1784
1785       subroutine schimc(p1x,p1y,p2x,p2y,out)
1786       implicit double precision(a-y)
1787       integer nx,jx,jy,nch,i1,i2,it,k
1788       double precision cc0(5),bm(5),bb0(5),pp0(5),bex(5),gaa(5)
1789       double precision sca(5,5,0:40000,2),sca1(5,5,0:40000,2)
1790       double precision a1(4),a2(4),q(4,20)
1791       complex*16 out,x0,x00
1792       common/pars/cc0,bm,bb0,pp0,bex,asp,sigo,gaa,ep,norm
1793       common/ipars/nch
1794       common/spac/sca,sca1
1795       common/mom/q
1796       common/wvars/sig0,bb,t1,t2
1797       common/alphas/alphap,alpha0,alphapr,alpha0r,alphapm,alpha0m
1798
1799       do k=1,4
1800          a1(k)=q(k,1)-q(k,3)
1801          a2(k)=q(k,2)-q(k,4)
1802       enddo
1803
1804       nx=10
1805       hx=2d0/dble(nx)
1806
1807       out=0d0
1808
1809       do jx=-nx,nx
1810          do jy=-nx,nx
1811
1812             tpx=dble(jx)*hx
1813             tpy=dble(jy)*hx
1814             tp2=tpx**2+tpy**2
1815             wt=hx*hx
1816
1817             p1xp=p1x-tpx
1818             p1yp=p1y-tpy
1819             t12=p1xp**2+p1yp**2
1820             p2xp=tpx+p2x
1821             p2yp=tpy+p2y
1822             t22=p2xp**2+p2yp**2
1823
1824             a1(1)=-p1xp
1825             a1(2)=-p1yp
1826             a2(1)=-p2xp
1827             a2(2)=-p2yp
1828
1829       do i1=1,nch
1830          do i2=1,nch
1831
1832            call screeningint(i1,i2,tp2,sc,sc1)
1833            call wev(x0,a1,a2)
1834
1835             x0=x0*pp0(i1)*pp0(i2)/dble(nch)**2
1836             x0=x0*dexp(-((t12+0.08d0+bb0(i1))*bex(i1))**cc0(i1)+
1837      &     (bex(i1)*(bb0(i1)+0.08d0))**cc0(i1))
1838             x0=x0*dexp(-((t22+0.08d0+bb0(i2))*bex(i2))**cc0(i2)+
1839      &     (bex(i2)*(bb0(i2)+0.08d0))**cc0(i2))
1840
1841             out=out+x0*wt*sc
1842
1843          enddo
1844       enddo
1845
1846
1847       enddo
1848       enddo
1849
1850       a1(1)=-p1x
1851       a1(2)=-p1y
1852       a2(1)=-p2x
1853       a2(2)=-p2y
1854
1855       t11=p1x**2+p1y**2
1856       t22=p2x**2+p2y**2
1857
1858       call formfac(t11,t22,x00p)
1859
1860       call wev(x00,a1,a2)
1861
1862
1863
1864       out=out+x00*x00p
1865
1866       return
1867       end
1868
1869 ccccc Pomeron -- diffrative eignstate i form factor
1870
1871       subroutine formfac(t1,t2,out)
1872       implicit double precision(a-y)
1873       double precision cc0(5),bm(5),bb0(5),pp0(5),bex(5),gaa(5)
1874       common/pars/cc0,bm,bb0,pp0,bex,asp,sigo,gaa,ep,norm
1875       integer nch,i1,i2
1876       common/ipars/nch
1877
1878       out=0d0
1879
1880       delta1=dexp(-t1/2d0)
1881       delta2=dexp(-t2/2d0)
1882
1883       do i1=1,nch
1884          do i2=1,nch
1885
1886             wt=gaa(i1)*gaa(i2)*pp0(i1)*pp0(i2)/dble(nch)**2
1887             wt=wt*dexp(-((t1+0.08d0+bb0(i1))*bex(i1))**cc0(i1)+
1888      &     (bex(i1)*(bb0(i1)+0.08d0))**cc0(i1))
1889             wt=wt*dexp(-((t2+0.08d0+bb0(i2))*bex(i2))**cc0(i2)+
1890      &     (bex(i2)*(bb0(i2)+0.08d0))**cc0(i2))
1891
1892             out=out+wt
1893
1894          enddo
1895       enddo
1896
1897       return
1898       end
1899
1900       subroutine screening(i,j,ktsq,out,out1)
1901       implicit double precision(a-y)
1902       integer i,j,ib,nb,it,nt,nch
1903       double precision cc0(5),bm(5),bb0(5),pp0(5),bex(5),gaa(5)
1904       common/pars/cc0,bm,bb0,pp0,bex,asp,sigo,gaa,ep,norm
1905       common/ipars/nch
1906       common/vars/s,rts,mmes,yx
1907
1908       pi=dacos(-1d0)
1909
1910       nb=5000
1911       hb=99d0/dble(nb)
1912
1913       out=0d0
1914       out1=0d0
1915
1916       do ib=1,nb
1917          bt=dble(ib)*hb
1918          wt=-bt/2d0/pi*hb
1919
1920          call opacityint(i,j,bt,fr,fr1)
1921
1922          sige=sigo*dexp(dlog(rts)*2d0*ep)
1923          fr=fr*gaa(i)*gaa(j)*sige
1924          fr1=fr1*gaa(i)*gaa(j)*sige
1925
1926          out=out+wt*(1d0-dexp(-fr/2d0))*besj0(bt*dsqrt(ktsq))
1927      &        *gaa(i)*gaa(j)
1928          out1=out1+wt*(1d0-dexp(-fr1/2d0))*besj0(bt*dsqrt(ktsq))
1929      &        *gaa(i)*gaa(j)
1930
1931       enddo
1932
1933       return
1934       end
1935
1936
1937       subroutine opacity(i,j,bt,out,out1)
1938       implicit double precision(a-y)
1939       integer i,j,ib,nb,it,nt,nch
1940       double precision cc0(5),bm(5),bb0(5),pp0(5),bex(5),gaa(5)
1941       common/pars/cc0,bm,bb0,pp0,bex,asp,sigo,gaa,ep,norm
1942       common/ipars/nch
1943       common/vars/s,rts,mmes,yx
1944
1945       pi=dacos(-1d0)
1946
1947       ampi=0.02d0
1948       amro=1d0
1949       a4=4d0*ampi
1950       alo=dlog(amro/ampi)
1951       bpol=2.4d0
1952
1953       nt=6000
1954       htt=6d0/dble(nt)**2
1955
1956          out=0d0
1957          out1=0d0
1958
1959          do it=0,nt
1960             t=dble(it)**2*htt
1961             if(it.eq.0) t=1d-8
1962             wt=htt*2d0*dble(it)/4d0/pi
1963             if(it.eq.0) wt=wt/2d0
1964             bes0=besj0(bt*dsqrt(t))
1965
1966             ffi=dexp(-((t+0.08d0+bb0(i))*bex(i))**cc0(i)+
1967      &           (bex(i)*(bb0(i)+0.08d0))**cc0(i))
1968             ffj=dexp(-((t+0.08d0+bb0(j))*bex(j))**cc0(j)+
1969      &           (bex(j)*(bb0(j)+0.08d0))**cc0(j))
1970
1971             asp1=asp
1972             form1=dlog(ffi*ffj)-2d0*t*asp1*dlog(rts)
1973 cccccccccccccccc
1974             ah=dsqrt(1d0+a4/t)
1975             h1pi=2d0*a4+t*(alo-ah*ah*ah*dlog((1d0+ah)/(ah-1d0))) !!! Pion loop insertion
1976             h1pi=h1pi*sigo/(72d0*pi**3*(1d0+t/bpol)**2)
1977 ccccccccccccccc
1978             ww=bes0*dexp(form1-2d0*h1pi*dlog(rts))
1979             aspt=t*asp+h1pi
1980
1981             out=out+ww*wt
1982             out1=out1+bes0*wt*ffi*ffj
1983
1984          enddo
1985
1986
1987       return
1988       end
1989
1990       function rgauss(r1,r2)
1991       implicit double precision (a-y)
1992
1993       pi=dacos(-1d0)
1994       rgauss=dsqrt(-2d0*dlog(r1))*dsin(2d0*pi*r2)
1995
1996       return
1997       end
1998
1999 c     boosting subroutine
2000
2001       subroutine boost(p,pboo,pcm,plb)
2002       real*8  pboo(4),pcm(4),plb(4),p,fact
2003          plb(4)=(pboo(4)*pcm(4)+pboo(3)*pcm(3)
2004      &             +pboo(2)*pcm(2)+pboo(1)*pcm(1))/p
2005          fact=(plb(4)+pcm(4))/(p+pboo(4))
2006          do 10 j=1,3
2007   10     plb(j)=pcm(j)+fact*pboo(j)
2008       return
2009       end
2010
2011 c     calculates lorentzian dot product for real 4-vectors
2012
2013       subroutine dot(v1,v2,dt)
2014       double precision v1(4),v2(4),dt
2015
2016       dt=-(v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3)
2017      &-v1(4)*v2(4))
2018
2019       return
2020       end
2021 c
2022
2023 *
2024 * subtractive mitchell-moore generator
2025 * ronald kleiss - october 2, 1987
2026 *
2027 * the algorithm is N(i)=[ N(i-24) - N(i-55) ]mod M,
2028 * implemented in a cirucular array with identifcation
2029 * of NR(i+55) and nr(i), such that effectively:
2030 *        N(1)   <---   N(32) - N(1)
2031 *        N(2)   <---   N(33) - N(2)  ....
2032 *   .... N(24)  <---   N(55) - N(24)
2033 *        N(25)  <---   N(1)  - N(25) ....
2034 *   .... N(54)  <---   N(30) - N(54)
2035 *        N(55)  <---   N(31) - N(55)
2036 *
2037 * in this version  M =2**30  and  RM=1/M=2.D0**(-30.D0)
2038 *
2039 * the array NR has been initialized by putting NR(i)=i
2040 * and subsequently running the algorithm 100,000 times.
2041 *
2042
2043       subroutine R2455(Ran)
2044       IMPLICIT REAL*8(A-H,O-Z)
2045       DIMENSION N(55)
2046       DATA N/
2047      . 980629335, 889272121, 422278310,1042669295, 531256381,
2048      . 335028099,  47160432, 788808135, 660624592, 793263632,
2049      . 998900570, 470796980, 327436767, 287473989, 119515078,
2050      . 575143087, 922274831,  21914605, 923291707, 753782759,
2051      . 254480986, 816423843, 931542684, 993691006, 343157264,
2052      . 272972469, 733687879, 468941742, 444207473, 896089285,
2053      . 629371118, 892845902, 163581912, 861580190,  85601059,
2054      . 899226806, 438711780, 921057966, 794646776, 417139730,
2055      . 343610085, 737162282,1024718389,  65196680, 954338580,
2056      . 642649958, 240238978, 722544540, 281483031,1024570269,
2057      . 602730138, 915220349, 651571385, 405259519, 145115737/
2058       DATA M/1073741824/
2059       DATA RM/0.9313225746154785D-09/
2060       DATA K/55/,L/31/
2061       IF(K.EQ.55) THEN
2062          K=1
2063       ELSE
2064          K=K+1
2065       ENDIF
2066       IF(L.EQ.55) THEN
2067          L=1
2068       ELSE
2069          L=L+1
2070       ENDIF
2071       J=N(L)-N(K)
2072       IF(J.LT.0) J=J+M
2073       N(K)=J
2074       RAN=J*RM
2075       END
2076
2077
2078
2079       SUBROUTINE RAMBO(N,ET,XM,P,WT)
2080 C------------------------------------------------------
2081 C
2082 C                       RAMBO
2083 C
2084 C    RA(NDOM)  M(OMENTA)  B(EAUTIFULLY)  O(RGANIZED)
2085 C
2086 C    A DEMOCRATIC MULTI-PARTICLE PHASE SPACE GENERATOR
2087 C    AUTHORS:  S.D. ELLIS,  R. KLEISS,  W.J. STIRLING
2088 C    THIS IS VERSION 1.0 -  WRITTEN BY R. KLEISS
2089 C
2090 C    N  = NUMBER OF PARTICLES (>1, IN THIS VERSION <101)
2091 C    ET = TOTAL CENTRE-OF-MASS ENERGY
2092 C    XM = PARTICLE MASSES ( DIM=100 )
2093 C    P  = PARTICLE MOMENTA ( DIM=(4,100) )
2094 C    WT = WEIGHT OF THE EVENT
2095 C
2096 C------------------------------------------------------
2097       implicit none
2098 !     boost arguments
2099       integer n
2100       double precision et,xm(100),p(4,100),wt
2101 !     function
2102       double precision rn
2103 !     internal variables
2104       double precision q(4,100),z(100),r(4),
2105      &     b(3),p2(100),xm2(100),e(100),v(100)
2106       integer iwarn(5)
2107       double precision acc
2108       integer itmax,ibegin
2109       data acc/1.0d-14/,itmax/6/,ibegin/0/,iwarn/5*0/
2110
2111       integer i,k,iter
2112       double precision twopi,po2log,xmt,c,s,f,rmas,g,a,x,bq,
2113      &     accu,f0,g0,wt2,wt3,wtm,x2,xmax
2114       integer nm
2115       save
2116
2117
2118 C
2119 C INITIALIZATION STEP: FACTORIALS FOR THE PHASE SPACE WEIGHT
2120       IF(IBEGIN.NE.0) GOTO 103
2121       IBEGIN=1
2122       TWOPI=8.*DATAN(1.D0)
2123       PO2LOG=DLOG(TWOPI/4.)
2124       Z(2)=PO2LOG
2125       DO 101 K=3,100
2126   101 Z(K)=Z(K-1)+PO2LOG-2.*DLOG(DFLOAT(K-2))
2127       DO 102 K=3,100
2128   102 Z(K)=(Z(K)-DLOG(DFLOAT(K-1)))
2129 C
2130 C CHECK ON THE NUMBER OF PARTICLES
2131   103 IF(N.GT.1.AND.N.LT.101) GOTO 104
2132       PRINT 1001,N
2133       STOP
2134 C
2135 C CHECK WHETHER TOTAL ENERGY IS SUFFICIENT; COUNT NONZERO MASSES
2136   104 XMT=0.
2137       NM=0
2138       DO 105 I=1,N
2139       IF(XM(I).NE.0.D0) NM=NM+1
2140   105 XMT=XMT+DABS(XM(I))
2141       IF(XMT.LE.ET) GOTO 201
2142       PRINT 1002,XMT,ET
2143       STOP
2144 C
2145 C THE PARAMETER VALUES ARE NOW ACCEPTED
2146 C
2147 C GENERATE N MASSLESS MOMENTA IN INFINITE PHASE SPACE
2148   201 DO 202 I=1,N
2149       C=2.*RN(1)-1.
2150       S=DSQRT(1.-C*C)
2151       F=TWOPI*RN(2)
2152       Q(4,I)=-DLOG(RN(3)*RN(4))
2153       Q(3,I)=Q(4,I)*C
2154       Q(2,I)=Q(4,I)*S*DCOS(F)
2155  202  Q(1,I)=Q(4,I)*S*DSIN(F)
2156
2157 C
2158 C CALCULATE THE PARAMETERS OF THE CONFORMAL TRANSFORMATION
2159       DO 203 I=1,4
2160   203 R(I)=0.
2161       DO 204 I=1,N
2162       DO 204 K=1,4
2163   204 R(K)=R(K)+Q(K,I)
2164       RMAS=DSQRT(R(4)**2-R(3)**2-R(2)**2-R(1)**2)
2165       DO 205 K=1,3
2166   205 B(K)=-R(K)/RMAS
2167       G=R(4)/RMAS
2168       A=1./(1.+G)
2169       X=ET/RMAS
2170 C
2171 C TRANSFORM THE Q'S CONFORMALLY INTO THE P'S
2172       DO 207 I=1,N
2173       BQ=B(1)*Q(1,I)+B(2)*Q(2,I)+B(3)*Q(3,I)
2174       DO 206 K=1,3
2175   206 P(K,I)=X*(Q(K,I)+B(K)*(Q(4,I)+A*BQ))
2176   207 P(4,I)=X*(G*Q(4,I)+BQ)
2177 C
2178 C CALCULATE WEIGHT AND POSSIBLE WARNINGS
2179       WT=PO2LOG
2180       IF(N.NE.2) WT=(2.*N-4.)*DLOG(ET)+Z(N)
2181       IF(WT.GE.-180.D0) GOTO 208
2182       IF(IWARN(1).LE.5) PRINT 1004,WT
2183       IWARN(1)=IWARN(1)+1
2184   208 IF(WT.LE. 174.D0) GOTO 209
2185       IF(IWARN(2).LE.5) PRINT 1005,WT
2186       IWARN(2)=IWARN(2)+1
2187 C
2188 C RETURN FOR WEIGHTED MASSLESS MOMENTA
2189   209 IF(NM.NE.0) GOTO 210
2190       WT=DEXP(WT)
2191       RETURN
2192 C
2193 C MASSIVE PARTICLES: RESCALE THE MOMENTA BY A FACTOR X
2194   210 XMAX=DSQRT(1.-(XMT/ET)**2)
2195       DO 301 I=1,N
2196       XM2(I)=XM(I)**2
2197   301 P2(I)=P(4,I)**2
2198       ITER=0
2199       X=XMAX
2200       ACCU=ET*ACC
2201   302 F0=-ET
2202       G0=0.
2203       X2=X*X
2204       DO 303 I=1,N
2205       E(I)=DSQRT(XM2(I)+X2*P2(I))
2206       F0=F0+E(I)
2207   303 G0=G0+P2(I)/E(I)
2208       IF(DABS(F0).LE.ACCU) GOTO 305
2209       ITER=ITER+1
2210       IF(ITER.LE.ITMAX) GOTO 304
2211 C      PRINT 1006,ITMAX,ACCU,DABS(F0)
2212       WRITE(99,1006)ITMAX,ACCU,DABS(F0)
2213       IF (DABS(F0).GT.1.0E-6) THEN
2214          WRITE(*,1007)DABS(F0)
2215       ENDIF
2216       GOTO 305
2217   304 X=X-F0/(X*G0)
2218       GOTO 302
2219   305 DO 307 I=1,N
2220       V(I)=X*P(4,I)
2221       DO 306 K=1,3
2222   306 P(K,I)=X*P(K,I)
2223   307 P(4,I)=E(I)
2224 C
2225 C CALCULATE THE MASS-EFFECT WEIGHT FACTOR
2226       WT2=1.
2227       WT3=0.
2228       DO 308 I=1,N
2229       WT2=WT2*V(I)/E(I)
2230   308 WT3=WT3+V(I)**2/E(I)
2231       WTM=(2.*N-3.)*DLOG(X)+DLOG(WT2/WT3*ET)
2232 C
2233 C RETURN FOR  WEIGHTED MASSIVE MOMENTA
2234       WT=WT+WTM
2235       IF(WT.GE.-180.D0) GOTO 309
2236       IF(IWARN(3).LE.5) PRINT 1004,WT
2237       IWARN(3)=IWARN(3)+1
2238   309 IF(WT.LE. 174.D0) GOTO 310
2239       IF(IWARN(4).LE.5) PRINT 1005,WT
2240       IWARN(4)=IWARN(4)+1
2241   310 WT=DEXP(WT)
2242       RETURN
2243 C
2244  1001 FORMAT(' RAMBO FAILS: # OF PARTICLES =',I5,' IS NOT ALLOWED')
2245  1002 FORMAT(' RAMBO FAILS: TOTAL MASS =',D15.6,' IS NOT',
2246      . ' SMALLER THAN TOTAL ENERGY =',D15.6)
2247  1004 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY UNDERFLOW')
2248  1005 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY  OVERFLOW')
2249  1006 FORMAT(' RAMBO WARNS:',I3,' ITERATIONS DID NOT GIVE THE',
2250      . ' DESIRED ACCURACY =',D10.3,' (STOPS AT',D10.3,')')
2251  1007 FORMAT(' RAMBO WARNS: END OF ITERATION ACCURACY TOO LOW =',D10.3)
2252       END
2253 C
2254
2255       function rn(idum)
2256 *
2257 * SUBTRACTIVE MITCHELL-MOORE GENERATOR
2258 * RONALD KLEISS - OCTOBER 2, 1987
2259 *
2260 * THE ALGORITHM IS N(I)=[ N(I-24) - N(I-55) ]MOD M,
2261 * IMPLEMENTED IN A CIRUCULAR ARRAY WITH IDENTIFCATION
2262 * OF NR(I+55) AND NR(I), SUCH THAT EFFECTIVELY:
2263 *        N(1)   <---   N(32) - N(1)
2264 *        N(2)   <---   N(33) - N(2)  ....
2265 *   .... N(24)  <---   N(55) - N(24)
2266 *        N(25)  <---   N(1)  - N(25) ....
2267 *   .... N(54)  <---   N(30) - N(54)
2268 *        N(55)  <---   N(31) - N(55)
2269 *
2270 * IN THIS VERSION  M =2**30  AND  RM=1/M=2.D0**(-30D0)
2271 *
2272 * THE ARRAY NR HAS BEEN INITIALIZED BY PUTTING NR(I)=I
2273 * AND SUBSEQUENTLY RUNNING THE ALGORITHM 100,000 TIMES.
2274 *
2275
2276       implicit none
2277       double precision rn
2278       integer idum
2279       integer n(55)
2280       data n/
2281      . 980629335, 889272121, 422278310,1042669295, 531256381,
2282      . 335028099,  47160432, 788808135, 660624592, 793263632,
2283      . 998900570, 470796980, 327436767, 287473989, 119515078,
2284      . 575143087, 922274831,  21914605, 923291707, 753782759,
2285      . 254480986, 816423843, 931542684, 993691006, 343157264,
2286      . 272972469, 733687879, 468941742, 444207473, 896089285,
2287      . 629371118, 892845902, 163581912, 861580190,  85601059,
2288      . 899226806, 438711780, 921057966, 794646776, 417139730,
2289      . 343610085, 737162282,1024718389,  65196680, 954338580,
2290      . 642649958, 240238978, 722544540, 281483031,1024570269,
2291      . 602730138, 915220349, 651571385, 405259519, 145115737/
2292       double precision eps
2293       double precision rm
2294       integer j,k,l,m
2295
2296       data eps/1D-9/
2297       data m/1073741824/
2298       data rm/0.9313225746154785D-09/
2299       data k/55/,l/31/
2300
2301
2302   1   CONTINUE
2303       IF(K.EQ.55) THEN
2304          K=1
2305       ELSE
2306          K=K+1
2307       ENDIF
2308       IF(L.EQ.55) THEN
2309          L=1
2310       ELSE
2311          L=L+1
2312       ENDIF
2313       J=N(L)-N(K)
2314       IF(J.LT.0) J=J+M
2315       N(K)=J
2316       RN=J*RM
2317       IF(RN.LT.EPS) GOTO 1
2318       IF(RN.GT.1D0-EPS) GOTO 1
2319       RETURN
2320       END