DIME/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,iin
  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, ntotal
  26 ccc new
  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
  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 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'
 129
 130 ccAM      rmax=1.8d1        ! max meson rapidity
 131 ccAM      rmin=-1.8d1       ! min meson rapidity
 132 ccAM      ecut=2d0          ! min meson p_t
 133
 134 c -- physics parameters
 135
 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"
 148           else
 149              write(6,*) "pp collisions"
 150           endif
 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,*) "*************************************************"
 164 cccccccc
 165
 166           ntotal=1000000             ! no. of runs for weighted events
 167           nev=100                    ! no. of unweighted events generated to event record
 168
 169 ccccccccc
 170
 171       if(formf.eq.'exp')then           ! Set parameters for Pomeron-Meson form factor in function 'fpi(x)'
 172          bexp=1d0/2.2d0
 173       elseif(formf.eq.'orexp')then
 174          bo=1d0/1.1d0
 175          ao=dsqrt(0.5d0)
 176       elseif(formf.eq.'power')then
 177          aoo=1.7d0
 178       endif
 179
 180 ccccccccccccccccccccccccccccccccccccccc
 181
 182 ccccccccc   Start of main body of code
 183
 184 cccccccccccccccccccccccccccccccccccccc
 185
 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
 191          sig0=13.63d0
 192          alphapm=0.7d0
 193          alpha0m=-0.7d0*mmes**2
 194          cf=31.79d0/0.389d0
 195          crho=4.23d0/0.389d0
 196
 197          nup=6
 198          istup(5)=1
 199          idup(5)=211
 200          mothup(1,5)=1
 201          mothup(2,5)=2
 202          icolup(1,5)=0
 203          icolup(2,5)=0
 204          vtimup(5)=0
 205          spinup(5)=9.
 206
 207          istup(6)=1
 208          idup(6)=-211
 209          mothup(1,6)=1
 210          mothup(2,6)=2
 211          icolup(1,6)=0
 212          icolup(2,6)=0
 213          vtimup(6)=0
 214          spinup(6)=9.
 215
 216       elseif(pflag.eq.'pi0')then
 217          mmes=0.1349766d0       ! pi0 mass, PDG 2011 value
 218          sig0=13.63d0
 219          alphapm=0.7d0
 220          alpha0m=-0.7d0*mmes**2
 221          cf=31.79d0/0.389d0
 222          crho=0d0
 223
 224          nup=6
 225          istup(5)=1
 226          idup(5)=111
 227          mothup(1,5)=1
 228          mothup(2,5)=2
 229          icolup(1,5)=0
 230          icolup(2,5)=0
 231          vtimup(5)=0
 232          spinup(5)=9.
 233
 234          istup(6)=1
 235          idup(6)=-111
 236          mothup(1,6)=1
 237          mothup(2,6)=2
 238          icolup(1,6)=0
 239          icolup(2,6)=0
 240          vtimup(6)=0
 241          spinup(6)=9.
 242
 243       elseif(pflag.eq.'kpkm')then
 244          mmes=0.493677d0         ! K+/- mass, PDG 2011 value
 245          sig0=11.82d0
 246          cf=17.255d0/0.389d0
 247          crho=9.105d0/0.389d0
 248
 249          nup=6
 250          istup(5)=1
 251          idup(5)=321
 252          mothup(1,5)=1
 253          mothup(2,5)=2
 254          icolup(1,5)=0
 255          icolup(2,5)=0
 256          vtimup(5)=0
 257          spinup(5)=9.
 258
 259          istup(6)=1
 260          idup(6)=-321
 261          mothup(1,6)=1
 262          mothup(2,6)=2
 263          icolup(1,6)=0
 264          icolup(2,6)=0
 265          vtimup(6)=0
 266          spinup(6)=9.
 267
 268       elseif(pflag.eq.'ks')then
 269          mmes=0.497614d0         ! K_0 mass, PDG 2011 value
 270          sig0=11.82d0
 271          cf=17.255d0/0.389d0
 272          crho=0d0
 273
 274          nup=6
 275          istup(5)=1
 276          idup(5)=310
 277          mothup(1,5)=1
 278          mothup(2,5)=2
 279          icolup(1,5)=0
 280          icolup(2,5)=0
 281          vtimup(5)=0
 282          spinup(5)=9.
 283
 284          istup(6)=1
 285          idup(6)=310
 286          mothup(1,6)=1
 287          mothup(2,6)=2
 288          icolup(1,6)=0
 289          icolup(2,6)=0
 290          vtimup(6)=0
 291          spinup(6)=9.
 292
 293       elseif(pflag.eq.'rho')then
 294          mmes0=0.77549d0          ! rho mass, PDG 2013 value
 295          mwidth=0.1491d0         ! rho width PDG 2013 value
 296          sig0=10d0
 297          cf=0d0
 298          crho=0d0
 299
 300          nup=10
 301          istup(5)=2
 302          idup(5)=113
 303          mothup(1,5)=1
 304          mothup(2,5)=2
 305          icolup(1,5)=0
 306          icolup(2,5)=0
 307          vtimup(5)=0
 308          spinup(5)=9.
 309
 310          istup(6)=2
 311          idup(6)=113
 312          mothup(1,6)=1
 313          mothup(2,6)=2
 314          icolup(1,6)=0
 315          icolup(2,6)=0
 316          vtimup(6)=0
 317          spinup(6)=9
 318
 319          do k=7,10
 320             istup(k)=1
 321             mothup(2,k)=0
 322             icolup(1,k)=0
 323             icolup(2,k)=0
 324             vtimup(k)=0
 325             spinup(k)=9.
 326          enddo
 327
 328          idup(7)=211
 329          idup(8)=-211
 330          idup(9)=211
 331          idup(10)=-211
 332          mothup(1,7)=5
 333          mothup(1,8)=5
 334          mothup(1,9)=6
 335          mothup(1,10)=6
 336
 337       endif
 338
 339 c -- other parameters
 340
 341       ebeam=rts/2d0
 342       s=rts*rts
 343       zi=(0d0,1d0)
 344       rt2=dsqrt(2d0)
 345       pi=dacos(-1d0)
 346       bp=rts/dsqrt(2.0d0)
 347       mp=0.93827d0
 348       beta=dsqrt(1d0-4d0*mp**2/s)
 349       s0=1d0
 350
 351 c     Pomeron + t-slope
 352
 353       bb=4d0
 354       bjac=6d0
 355       bjac1=2d0
 356
 357       alphap=0.25d0    ! D-L '92 fit
 358       alpha0=1.0808d0
 359       alphapr=0.93d0
 360       alpha0r=0.55d0
 361
 362       mf127=1.275d0
 363       mf1525=1.525d0
 364
 365       cpom=sig0/0.389d0
 366       aff=-0.860895d0
 367       ar=-1.16158d0
 368
 369       mmes1=mmes
 370       mmes2=mmes
 371
 372 ccccc Initialise rambo (rho0 decay)
 373
 374       if(pflag.eq.'rho')then
 375
 376          mmes=mmes0
 377          npart=2
 378          do j=1,4
 379             am(j)=0.13957018d0
 380          enddo
 381       endif
 382
 383 ccccccc
 384
 385 c     initialise counters etc
 386
 387 c      do hh=1,20
 388
 389       nhist=1
 390       sum=0d0
 391       sum1=0d0
 392       ncut=0
 393
 394       weightm=0d0
 395
 396 c     initialise histograms
 397
 398       do i=1,nhist
 399          call histo3(i)
 400       enddo
 401
 402       num=0
 403
 404 cccccccccccccccccccccccccccccccccccccccccccccccccccccc
 405
 406 c     Incoming protons to event record array
 407
 408 cccccccccccccccccccccccccccccccccccccccccccccccccccccc
 409
 410       ID(1)=2212
 411       q(1,1)=0d0
 412       q(2,1)=0d0
 413       q(3,1)=ebeam*beta
 414       q(4,1)=ebeam
 415
 416       istup(1)=-1
 417       idup(1)=2212
 418       mothup(1,1)=0
 419       mothup(2,1)=0
 420       icolup(1,1)=0
 421       icolup(2,1)=0
 422       do j=1,4
 423          pup(j,1)=q(j,1)
 424       enddo
 425       pup(5,1)=dsqrt(q(4,1)**2-q(3,1)**2-q(2,1)**2-q(1,1)**2)
 426       vtimup(1)=0
 427       spinup(1)=9.
 428
 429       q(1,2)=0d0
 430       q(2,2)=0d0
 431       q(3,2)=-ebeam*beta
 432       q(4,2)=ebeam
 433
 434       istup(2)=-1
 435       if(ppbar.eq.'true')then
 436          idup(2)=-2212
 437       else
 438          idup(2)=2212
 439       endif
 440       mothup(1,2)=0
 441       mothup(2,2)=0
 442       icolup(1,2)=0
 443       icolup(2,2)=0
 444       do j=1,4
 445          pup(j,2)=q(j,2)
 446       enddo
 447       pup(5,2)=dsqrt(q(4,2)**2-q(3,2)**2-q(2,2)**2-q(1,2)**2)
 448       vtimup(2)=0
 449       spinup(2)=9.
 450
 451 cccc  outgoing initial info
 452
 453       istup(3)=1
 454       if(ppbar.eq.'true')then
 455          idup(2)=-2212
 456       else
 457          idup(2)=2212
 458       endif
 459       idup(3)=2212
 460       mothup(1,3)=1
 461       mothup(2,3)=2
 462       icolup(1,3)=0
 463       icolup(2,3)=0
 464       vtimup(3)=0
 465       spinup(3)=9
 466
 467       istup(4)=1
 468       idup(4)=2212
 469       mothup(1,4)=1
 470       mothup(2,4)=2
 471       icolup(1,4)=0
 472       icolup(2,4)=0
 473       vtimup(4)=0
 474       spinup(4)=9
 475
 476 ccc   HEPEVT
 477       if(output.eq.'hepevt')then
 478
 479       nhep=nup
 480       do k=1,5
 481          phep(k,1)=pup(k,1)
 482          phep(k,2)=pup(k,2)
 483       enddo
 484
 485       do k=1,nup
 486          isthep(k)=istup(k)
 487          idhep(k)=idup(k)
 488          jmohep(1,k)=mothup(1,k)
 489          jmohep(2,k)=jmohep(1,k)
 490          jdahep(1,k)=0
 491          jdahep(2,k)=0
 492          vhep(1,k)=0d0
 493          vhep(2,k)=0d0
 494          vhep(3,k)=0d0
 495          vhep(4,k)=0d0
 496       enddo
 497
 498       if(pflag.eq.'rho')then
 499       jdahep(1,5)=7
 500       jdahep(2,5)=8
 501       jdahep(1,6)=9
 502       jdahep(2,6)=10
 503       endif
 504
 505       jmohep(2,5)=2
 506       jmohep(2,6)=2
 507
 508       endif
 509
 510       if(unw.eq.'true')then
 511          lmax=2
 512       else
 513          lmax=1
 514       endif
 515 c
 516 c
 517 c     first  generate weighted events first
 518       ll = 1
 519       ip=ntotal+1
 520       write(6,*)'Generating weighted events...'
 521       do i=1, ntotal
 522          isuc = 0
 523          call dimegenerate(isuc)
 524       enddo
 525 c     prepare for unweighted generation
 526       ll = 2
 527       write(6,*)'Generating unweighted events...'
 528 c
 529       return
 530       end
 531
 532       subroutine dimegenerate(isuccess)
 533       implicit none
 534
 535 ccccc hepevt output
 536       integer nmxhep,kk
 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
 543       INTEGER MAXNUP
 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),
 550      &              SPINUP(MAXNUP)
 551 ccccc rambo variables
 552       integer npart,nrun
 553       double precision pin(4), am(100), pout(4,100), wt, ein
 554       common/ramboc/ pin, am, pout, wt, ein,
 555      &npart, nrun
 556
 557 ccccc
 558       double precision etaelmax,etaelmin,ptelmin,ptphmin,ecut,rmax,
 559      &rmin,mcut
 560       common/cuts/etaelmax,etaelmin,ptelmin,ptphmin,ecut,rmax,rmin,mcut
 561 ccccc
 562       double precision q(4,20)
 563       common /mom/ q
 564 ccccc
 565       double precision pt1(2), pt2(2), ptx(2)
 566       common/mompt/pt1,pt2,ptx
 567 ccccc
 568       double precision s, rts, mmes, yx
 569       double precision bjac, bp
 570       double precision mf127, mf1525
 571
 572       double precision m0, mmes0, mmes1, mmes2, mp, mwidth, pi, rt2,
 573      &ebeam, sh, th, uh, sum, sum1, weightm
 574
 575       integer id(20), ncut, ll, icut, nev, num, iin, ntotal
 576
 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,
 581      &              bjac, bp
 582       common /ivars/ ncut, ll, icut, nev, num, ntotal
 583 ccccc
 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
 587       common/ff/formf
 588 ccccc
 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,
 602      & yx2
 603       double precision pboo(4), pcm(4), plb(4)
 604       double precision svec(4)
 605       complex*16 zmat
 606       integer h, i, j, k, isuccess
 607       weight=0d0
 608       isuccess = 0
 609       call r2455(ran0)
 610       call r2455(ran1)
 611       call r2455(ran2)
 612       call r2455(ran3)
 613       call r2455(ran4)
 614       call r2455(ran5)
 615       call r2455(ran6)
 616
 617 ccccccccccccccccccccccccccccccccccccccccccccccccccc
 618
 619 c     incoming protons
 620
 621       ID(1)=2212
 622       q(1,1)=0d0
 623       q(2,1)=0d0
 624       q(3,1)=ebeam
 625       q(4,1)=ebeam
 626
 627       ID(2)=2212
 628       q(1,2)=0d0
 629       q(2,2)=0d0
 630       q(3,2)=-ebeam
 631       q(4,2)=ebeam
 632
 633       phi1=2d0*pi*ran0
 634       phi2=2d0*pi*ran1
 635
 636       pt1sq=-dlog(ran2)/(bjac)
 637       pt2sq=-dlog(ran3)/(bjac)
 638
 639       weight=(dexp(bjac*pt1sq)*dexp(bjac*pt2sq))/bjac**2
 640
 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)
 645
 646       pt1x=pt1(1)
 647       pt1y=pt1(2)
 648       pt2x=pt2(1)
 649       pt2y=pt2(2)
 650
 651 cccccccccccccccccccccccccccccccc
 652
 653 ccc non-zero rho width
 654
 655 cccccccccccccccccccccccccccccccc
 656
 657        if(pflag.eq.'rho')then  ! new part here
 658
 659  677     call r2455(rm1)
 660          call r2455(rm2)
 661
 662          msmax=mmes0+4d0*mwidth
 663          msmin=2d0*0.13957018d0
 664
 665 cccccccccccccccc
 666
 667          almin=datan(-(mmes0**2-msmin**2)/mwidth/mmes0)
 668          almax=datan(-(mmes0**2-msmax**2)/mwidth/mmes0)
 669
 670          al1=almin+(almax-almin)*rm1
 671          al2=almin+(almax-almin)*rm2
 672
 673          mmes1=dsqrt(dtan(al1)*mmes0*mwidth+mmes0**2)
 674          mmes2=dsqrt(dtan(al2)*mmes0*mwidth+mmes0**2)
 675
 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
 680
 681 ccccccccccccccccc
 682
 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
 687
 688          if(mmes1.lt.2d0*0.13957018d0) goto 677
 689          if(mmes2.lt.2d0*0.13957018d0) goto 677
 690
 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)
 698
 699       endif
 700
 701
 702 cccccccccccccccccccccccccccccccccccccccccccccccccccc
 703
 704 c     pi rapidities
 705
 706          call r2455(ranx1)
 707          call r2455(ranx2)
 708          call r2455(ranx3)
 709          call r2455(ranx4)
 710
 711          phix1=2d0*pi*ranx1
 712
 713 c         r2max=1d0
 714 c         r2=r2max*ranx1
 715 c         ptxsq1=-dlog(r2)/bjac1  ! generate exp p_t^2 (can be more efficient)
 716
 717          ptxsqmin=0d0
 718          ptxsqmax=10d0         ! generate linear p_t^2
 719
 720          ptxsq1=ptxsqmin+(ptxsqmax-ptxsqmin)*ranx2
 721
 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
 727
 728          rmx1=dsqrt(mmes1**2+ptxsq1)
 729          rmx2=dsqrt(mmes2**2+ptxsq2)
 730
 731          ymax1=dlog(rts/rmx1)
 732
 733          if(cuts.eq.'true')then
 734             if(rmax.lt.ymax1)then
 735                ymax1=rmax
 736             endif
 737          endif
 738
 739          ymin1=-ymax1
 740
 741          yx1=ymin1+(ymax1-ymin1)*ranx3
 742
 743          ymax2=(dlog((rts-rmx1*dexp(yx1))/rmx2))
 744          ymin2=-(dlog((rts-rmx1*dexp(-yx1))/rmx2))
 745
 746          if(cuts.eq.'true')then
 747             if(ymax2.gt.rmax)then
 748                ymax2=rmax
 749             endif
 750             if(ymin2.lt.-rmax)then
 751                ymin2=-rmax
 752             endif
 753          endif
 754
 755          if(ymax2.lt.ymin2)then
 756             weight=0d0
 757             goto 700
 758          endif
 759
 760          yx2=ymin2+(ymax2-ymin2)*ranx4
 761          x1=(rmx2*dexp(yx2)+rmx1*dexp(yx1))/rts
 762          x2=(rmx2*dexp(-yx2)+rmx1*dexp(-yx1))/rts
 763
 764          weight=weight*(ptxsqmax-ptxsqmin)
 765 c         weight=weight*dexp(bjac1*ptxsq1)*r2max/bjac1
 766          weight=weight*(ymax1-ymin1)*(ymax2-ymin2)
 767
 768
 769          q(3,5)=rmx1*(dexp(yx1)-dexp(-yx1))/2d0
 770          q(4,5)=rmx1*(dexp(yx1)+dexp(-yx1))/2d0
 771
 772          q(3,6)=rmx2*(dexp(yx2)-dexp(-yx2))/2d0
 773          q(4,6)=rmx2*(dexp(yx2)+dexp(-yx2))/2d0
 774
 775 ccccccccccccccccccccccccccccccccccccccccccccccccccccc
 776
 777 c     impose massive on-shell condition by solving
 778 c                   p1+ + cc1/p2- = aa1
 779 c                   p2- + cc2/p1+ = aa2
 780
 781       aa1=bp*(1d0-x1)
 782       aa2=bp*(1d0-x2)
 783       cc1=0.5d0*(pt2sq+mp**2)
 784       cc2=0.5d0*(pt1sq+mp**2)
 785
 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
 789          weight=0d0
 790          goto 700
 791       endif
 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)
 796
 797       if(p1p.lt.0d0.or.p1m.lt.0d0.or.p2p.lt.0d0.or.p2m.lt.0d0)then
 798          weight=0d0
 799          goto 700
 800       endif
 801
 802       t1=-rts*p1m*rt2
 803       t2=-rts*p2p*rt2
 804
 805       q(1,3)=pt1(1)
 806       q(2,3)=pt1(2)
 807       q(3,3)=(p1p-p1m)/rt2
 808       q(4,3)=(p1p+p1m)/rt2
 809
 810       q(1,4)=pt2(1)
 811       q(2,4)=pt2(2)
 812       q(3,4)=(p2p-p2m)/rt2
 813       q(4,4)=(p2p+p2m)/rt2
 814
 815 ccccccccccccccccccccccccccccccccccccccccccccccccc
 816
 817       do k=1,4
 818          svec(k)=q(k,5)+q(k,6)
 819       enddo
 820
 821          call dot(svec,svec,sh)
 822
 823 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 824
 825 c     rho0 decays
 826
 827 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 828
 829          if(pflag.eq.'rho')then
 830
 831             kk=6
 832
 833             do k=5,6
 834
 835                if(k.eq.5)then
 836                   mmesp=mmes1
 837                else
 838                   mmesp=mmes2
 839                endif
 840                call rambo(npart,mmesp,am,pout,wt) ! call RAMBO
 841
 842                do j=1,4
 843                   pboo(j)=q(j,k)
 844                enddo
 845
 846                do h=1,2
 847                   do j=1,4
 848                      pcm(j)=pout(j,h)
 849                   enddo
 850                   call boost(mmesp,pboo,pcm,plb)
 851                   kk=kk+1
 852                   do j=1,4
 853                      q(j,kk)=plb(j)
 854                   enddo
 855                enddo
 856
 857             enddo
 858
 859          endif
 860
 861
 862 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 863
 864 c      place cuts
 865
 866          if(cuts.eq.'true')then
 867             call cut(icut)
 868
 869             if(icut.eq.0)then
 870                weight=0d0
 871                weight0=0d0
 872                weight1=0d0
 873                weight2=0d0
 874                ncut=ncut+1
 875                goto 700
 876             endif
 877
 878          endif
 879
 880 c         print*,i
 881
 882 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 883
 884 c     Write 4-momenta to event record array
 885
 886 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 887
 888          if(output.eq.'lhe')then
 889
 890          do k=3,4
 891             do j=1,4
 892                pup(j,k)=q(j,k)
 893             enddo
 894             pup(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2)
 895          enddo
 896
 897          if(pflag.eq.'rho')then
 898
 899             do k=5,10
 900                do j=1,4
 901                   pup(j,k)=q(j,k)
 902                enddo
 903                pup(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2)
 904             enddo
 905
 906          else
 907
 908             do k=5,6
 909                do j=1,4
 910                   pup(j,k)=q(j,k)
 911                enddo
 912                pup(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2)
 913             enddo
 914
 915          endif
 916
 917          else
 918
 919          do k=3,4
 920             do j=1,4
 921                phep(j,k)=q(j,k)
 922             enddo
 923             phep(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2)
 924          enddo
 925
 926          if(pflag.eq.'rho')then
 927
 928             do k=5,10
 929                do j=1,4
 930                   phep(j,k)=q(j,k)
 931                enddo
 932                phep(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2)
 933             enddo
 934
 935          else
 936
 937             do k=5,6
 938                do j=1,4
 939                   phep(j,k)=q(j,k)
 940                enddo
 941                phep(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2)
 942             enddo
 943
 944          endif
 945
 946          endif
 947
 948 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 949
 950 c     Weight evaluation
 951
 952 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 953
 954 c     pp-->pp(M_1 M_2) matrix element
 955
 956 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 957
 958
 959       call schimc(pt1x,pt1y,pt2x,pt2y,zmat)
 960
 961       weight=weight*cdabs(zmat)**2
 962       weight=weight/norm**2
 963
 964 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 965
 966       weight=weight/(s**2*(16d0)**3*pi**5)*389d3
 967
 968 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 969
 970 ccc   Probability of no additional particles produced in production subprocess
 971
 972 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 973
 974       if(fsi.eq.'true')then
 975
 976          if(pflag.eq.'pipm'.or.pflag.eq.'pi0')then
 977
 978             if(dsqrt(sh).gt.mf127)then
 979                m0=mf127
 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          elseif(pflag.eq.'kpkm'.or.pflag.eq.'ks')then
 989
 990             if(dsqrt(sh).gt.mf1525)then
 991                m0=mf1525
 992                c=0.7d0
 993                cont=1d0-(m0-2d0*mmes)**2/m0**2
 994                exn=c*dlog(((dsqrt(sh)-2d0*mmes)**2)/m0**2+cont)
 995                snp=dexp(-exn)
 996             else
 997                snp=1d0
 998             endif
 999
1000          else
1001
1002             snp=1d0
1003
1004          endif
1005
1006       else
1007
1008          snp=1d0
1009
1010       endif
1011
1012       if(snp.gt.1d0)then
1013          snp=1d0
1014       endif
1015
1016       weight=weight*snp
1017
1018 ccccccccccccccc
1019
1020       if(pflag.eq.'pi0'.or.pflag.eq.'ks'.or.pflag.eq.'rho')then
1021          weight=weight/2d0    ! symmetry factor
1022       endif
1023
1024 cccccccccccccccccccccccccccccccccccccccccccccccccccccc
1025
1026  666  weight=weight
1027
1028       if(weightm.lt.weight)then
1029          weightm=weight
1030       endif
1031
1032
1033       if(ll.eq.1)then
1034
1035          xwgtup=weight
1036
1037          wthist=weight/dfloat(ntotal)
1038
1039          call binit(wthist)
1040
1041       else
1042
1043          call r2455(ranhist)
1044
1045          if(ranhist.lt.weight/weightm)then
1046             isuccess = 1
1047             xwgtup=1d0
1048
1049             num=num+1
1050
1051
1052 c            call binit(sumt/nev)
1053
1054             write(35,302) num, nup
1055
1056             if(output.eq.'lhe')then
1057
1058             do j=1,nup
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)
1062             enddo
1063
1064             else
1065
1066                do j=1,nup
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)
1071                enddo
1072
1073             endif
1074
1075
1076          endif
1077
1078       endif
1079
1080       if(ll.eq.2)then
1081
1082       if(num.gt.nev-1)then  ! exit loop once required no. of unweighted events generated
1083
1084          ntotal=i
1085          return
1086       endif
1087       endif
1088
1089  700  sum=sum+weight
1090       sum1=sum1+weight**2
1091
1092 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1093 c                                                                    c
1094 c     End of event loop                                              c
1095 c                                                                    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
1102      &,E13.6)
1103  302  format(1x,i8,1x,i5,1x,E13.6)
1104
1105       return
1106       end
1107
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
1111
1112  888  continue
1113       sum=sum/dfloat(ntotal)
1114       sum1=sum1/dfloat(ntotal)
1115       var=dsqrt((sum1-sum**2)/dfloat(ntotal))
1116       EFF=1d0*(ntotal-ncut)/ntotal
1117
1118       if(ll.eq.1)then
1119          print*,'...done!'
1120          write(6,*)
1121          sumt=sum
1122       else
1123          print*,'...done!'
1124          write(6,*)
1125       endif
1126
1127
1128
1129 ccccc     create histograms
1130
1131       if(ll.eq.1)then
1132
1133       do i=1,nhist
1134            call histo2(i,0)
1135       enddo
1136
1137       endif
1138
1139       if(ll.eq.1)then
1140       if(pflag.eq.'pipm')then
1141          write(6,*)
1142          write(6,*)'pi+pi- production'
1143       elseif(pflag.eq.'pi0')then
1144          write(6,*)
1145          write(6,*)'pi0pi0 production'
1146       elseif(pflag.eq.'kpkm')then
1147          write(6,*)
1148          write(6,*)'K+K- production'
1149       elseif(pflag.eq.'ks')then
1150          write(6,*)
1151          write(6,*)'K0K0 production'
1152       elseif(pflag.eq.'rho')then
1153          write(6,*)
1154          write(6,*)'rho0rho0 production'
1155       endif
1156       endif
1157
1158       if(ll.eq.1)then
1159          if(pflag.eq.'rho')then
1160          write(6,221)rts,ntotal,sum,var
1161       else
1162          write(6,222)rts,mmes,ntotal,sum,var
1163       endif
1164       else
1165          if(output.eq.'lhe')then
1166             print*,'LHE output'
1167          else
1168             print*,'HEPEVT output'
1169          endif
1170          write(6,223)nev,sum,var
1171       endif
1172
1173 !      enddo
1174
1175       close(35)
1176
1177  221  format(/,
1178      .       3x,'collider energy  (GeV)    : ',f10.3,/,
1179      .       3x,'number of events          : ',i9,/,
1180      .       3x,'sigma (nb)                : ',f15.5,'  +-  ',f15.5,
1181      .   /)
1182
1183  222  format(/,
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,
1188      .   /)
1189
1190  223  format(3x,'number of events          : ',i9,/,
1191      .       3x,'sigma (nb)                : ',f15.5,'  +-  ',f15.5,
1192      .   /)
1193
1194       call cpu_time(t2)
1195       print*,'time elapsed = ', t2, ' s'
1196
1197       return
1198       end
1199
1200 ccccc Pomeron -- (off-shell) meson form factor
1201
1202       function fpi(x)
1203       implicit double precision(a-y)
1204       character formf*10
1205       common/vars/s,rts,mmes,yx, iin
1206       common/ff/formf
1207       common/ffpars/bexp,ao,bo,aoo,boo
1208
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)
1215       endif
1216
1217       return
1218       end
1219
1220 ccccccc  Pom Pom --> meson pair amplitude
1221
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)
1226       integer k, iin
1227       character*10 mregge,pflag
1228       complex*16 matt
1229       common/mom/q
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
1235       common/regge/mregge
1236       common/flags/ pflag, fsi, ppbar, output, cuts, unw
1237
1238
1239       zi=(0d0,1d0)
1240       pi=dacos(-1d0)
1241
1242       do k=1,4
1243          q(k,10)=v1(k)
1244          q(k,11)=v2(k)
1245       enddo
1246
1247       do k=1,4
1248          tvec(k)=q(k,10)-q(k,5)
1249          uvec(k)=q(k,10)-q(k,6)
1250       enddo
1251
1252          call dot(tvec,tvec,th)
1253          call dot(uvec,uvec,uh)
1254
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
1257
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
1260
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
1263
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
1266
1267        bt1=(2d0*alphap*dlog(sht1))
1268        bt2=(2d0*alphap*dlog(sht2))
1269        bu1=(2d0*alphap*dlog(shu1))
1270        bu2=(2d0*alphap*dlog(shu2))
1271
1272        bt1r=(2d0*alphapr*dlog(sht1))
1273        bt2r=(2d0*alphapr*dlog(sht2))
1274        bu1r=(2d0*alphapr*dlog(shu1))
1275        bu2r=(2d0*alphapr*dlog(shu2))
1276
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
1281
1282        t11=v1(1)**2+v1(2)**2
1283        t22=v2(1)**2+v2(2)**2
1284
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)
1292
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))
1296
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)
1301
1302           matt=zmu+zmt
1303        return
1304        end
1305
1306 c     binning subroutine
1307
1308       subroutine binit(wt)
1309       implicit double precision(a-y)
1310       double precision q(4,20),pt1(2),pt2(2),ptx(2)
1311       common/mom/q
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
1318
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)))
1321
1322       ypi1=0.5d0*dlog((q(4,5)+q(3,5))
1323      &     /(q(4,5)-q(3,5)))
1324       ypi2=0.5d0*dlog((q(4,6)+q(3,6))
1325      &     /(q(4,6)-q(3,6)))
1326
1327       call histo1(1,30,0d0,4.5d0,dsqrt(sh),wt)
1328
1329       return
1330       end
1331
1332 ccccc Cutting subroutine
1333 ccccc
1334 ccccc User can define their own cuts on any particle momenta, stored in array q(i,j)
1335 ccccc
1336
1337       subroutine cut(icut)
1338       implicit double precision(a-y)
1339       double precision q(4,20)
1340       integer icut
1341       common/mom/q
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
1347       integer iii
1348       common/ivar/iii
1349
1350       icut=0
1351
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)))
1357
1358       ypi1=.5d0*dlog((q(4,5)+q(3,5))
1359      .              /(q(4,5)-q(3,5)))
1360
1361         etpi1=q(4,5)*dsqrt(pmodpi1**2-q(3,5)**2)/pmodpi1
1362
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)))
1368
1369       ypi2=.5d0*dlog((q(4,6)+q(3,6))
1370      .              /(q(4,6)-q(3,6)))
1371
1372         etpi2=q(4,6)*dsqrt(pmodpi2**2-q(3,6)**2)/pmodpi2
1373
1374 c        print*,ptpi1,ptpi2
1375
1376 ccccc    example cuts....
1377
1378 c        if(dsqrt(sh).lt.0.8d0)return
1379
1380         xf1=2d0*q(3,3)/rts
1381         xf2=2d0*q(3,4)/rts
1382
1383 c        if(dabs(xf1).lt.0.9d0)return
1384 c        if(dabs(xf2).lt.0.9d0)return
1385
1386 ccccccccc  Basic cuts described at start of code
1387
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
1394
1395         icut=1
1396
1397       return
1398       end
1399
1400 c     prints histograms
1401
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 /'*'/
1408       save
1409       open(10,file="output.dat")
1410       y0(ih)=x0
1411       y1(ih)=x1
1412       ic(ih)=ib
1413       if(x.lt.x0) goto 11
1414       if(x.gt.x1) goto 12
1415       ix=idint((x-x0)/(x1-x0)*dble(ib))+1
1416       h(ih,ix)=h(ih,ix)+w
1417       if(h(ih,ix).gt.hx(ih)) hx(ih)=h(ih,ix)
1418       ii(ih)=ii(ih)+1
1419       return
1420    11 iu(ih)=iu(ih)+1
1421       return
1422    12 io(ih)=io(ih)+1
1423       return
1424       entry histo2(ih,il)
1425       ib1=ic(ih)
1426       x01=y0(ih)
1427       x11=y1(ih)
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
1435       write(6,23)
1436    23 format(35(1h ),3(10h----+----i))
1437       do 27 iv=1,ib1
1438       z=(dble(iv)-0.5d0)/dble(ib1)*(x11-x01)+x01
1439       if(il.eq.1) goto 24
1440       iz=idint(h(ih,iv)/hx(ih)*30.)+1
1441       goto 25
1442    24 iz=-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
1453    27 continue
1454       write(6,23)
1455       return
1456    28 write(6,29)
1457    29 format('  no entries inside histogram')
1458       return
1459       entry histo3(ih)
1460       do 31 i=1,100
1461    31 h(ih,i)=0.
1462       hx(ih)=0.
1463       io(ih)=0
1464       iu(ih)=0
1465       ii(ih)=0
1466       close(10)
1467       return
1468       end
1469
1470
1471 cccc  Initializes soft model parameters
1472
1473       subroutine initpars(in)
1474       implicit double precision(a-y)
1475       integer in,i1,i2
1476       integer nch
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
1479       common/ipars/nch
1480       common/vars/s,rts,mmes,yx, iin
1481
1482       pi=dacos(-1d0)
1483       if(in.eq.1)then
1484
1485          ep=0.13d0
1486          asp=0.08d0
1487          ep1=0d0
1488          sigo=23d0/0.39d0
1489          gaa(3)=0.4d0
1490          gd2=0.3d0
1491          nch=2
1492          ntf=0
1493
1494          cc0(1)=0.45d0
1495          bm(1)=3d0
1496          bb0(1)=0.1d0
1497          pp0(1)=0.92d0
1498          bex(1)=8.5d0
1499
1500          cc0(2)=0.45d0
1501          bm(2)=1.5d0
1502          bb0(2)=0.5d0
1503          pp0(2)=0.1d0
1504          bex(2)=4.5d0
1505
1506          cc0(3)=1d0
1507          bm(3)=0d0
1508          bb0(3)=0.8d0
1509          pp0(3)=0.5d0
1510          bex(3)=0.5d0
1511
1512       elseif(in.eq.2)then
1513
1514          ep=0.115d0
1515          asp=0.11d0
1516          ep1=0d0
1517          sigo=33d0/0.39d0
1518          gaa(3)=0.6d0
1519          gd2=0.16d0
1520          nch=2d0
1521          ntf=0d0
1522
1523          cc0(1)=0.63d0
1524          bm(1)=3d0
1525          bb0(1)=0.1d0
1526          pp0(1)=0.5d0
1527          bex(1)=8d0
1528
1529          cc0(2)=0.47d0
1530          bm(2)=1.5d0
1531          bb0(2)=0.5d0
1532          pp0(2)=0.1d0
1533          bex(2)=6d0
1534
1535          cc0(3)=1d0
1536          bm(3)=0d0
1537          bb0(3)=0.8d0
1538          pp0(3)=0.5d0
1539          bex(3)=0.5d0
1540
1541       elseif(in.eq.3)then
1542
1543          ep=0.093d0
1544          asp=0.075d0
1545          ep1=0d0
1546          sigo=60d0/0.39d0
1547          gaa3=4.8d0
1548          gd2=1.03d0
1549          nch=2
1550          nga=1
1551          cc0(1)=0.55d0
1552          bm(1)=3d0
1553          bb0(1)=0.27d0
1554          pp0(1)=0.48d0
1555          bex(1)=5.3d0
1556
1557          cc0(2)=0.48d0
1558          bm(2)=1.5d0
1559          bb0(2)=0.1d0
1560          pp0(2)=1d0
1561          bex(2)=3.8d0
1562
1563          cc0(3)=0.24d0
1564
1565       elseif(in.eq.4)then
1566
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)
1573          nch=2
1574          nga=1                  !!! A flag (?)
1575          cc0(1)=0.6d0           !!! d1
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
1579          bex(1)=7.2d0           !!! b1
1580
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
1585          bex(2)=4.2d0           !!! b2
1586
1587          cc0(3)=0.12d0          !!! Beta : k^2_min ~ s^Beta
1588
1589       endif
1590
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
1594
1595       if(in.eq.3.or.in.eq.4)then
1596
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)
1602
1603       elseif(in.eq.1.or.in.eq.2)then
1604
1605          if(nch.eq.2)then
1606             gaa(1)=1d0+dsqrt(gd2)
1607             gaa(2)=1d0-dsqrt(gd2)
1608             gaa(3)=0.
1609          elseif(nch.eq.1)then
1610             gaa(1)=1d0
1611             gaa(2)=0d0
1612             gaa(3)=0d0
1613          endif
1614
1615       endif
1616
1617       sum=0d0
1618
1619       do i1=1,nch
1620          do i2=1,nch
1621             sum=sum+gaa(i1)*gaa(i2)*pp0(i1)*pp0(i2)/dble(nch)**2
1622          enddo
1623       enddo
1624
1625       norm=sum
1626       return
1627       end
1628
1629 ccccc Calculates proton opacity
1630
1631       subroutine calcop
1632       implicit double precision(a-y)
1633       integer nb,i1,i2,nch,ib
1634       common/ipars/nch
1635       double precision op(5,5,10000,2),oph(5,5,10000,2)
1636       common/opac/op,oph
1637
1638       nb=900
1639       hb=100d0/dble(nb)
1640
1641       print*,'Calculating opacity'
1642
1643
1644       do ib=1,nb+1
1645          bt=dble(ib-1)*hb
1646
1647          do i1=1,nch
1648             do i2=1,nch
1649
1650                call opacity(i1,i2,bt,fr,fr1)
1651
1652
1653                op(i1,i2,ib,1)=bt
1654                op(i1,i2,ib,2)=fr
1655                oph(i1,i2,ib,1)=bt
1656                oph(i1,i2,ib,2)=fr1
1657
1658                write(40,*)bt,fr,fr1
1659
1660             enddo
1661          enddo
1662       enddo
1663
1664       return
1665       end
1666
1667 cccc  Calculates screening amplitude (in k_t space)
1668
1669       subroutine calcscreen
1670       implicit double precision(a-y)
1671       integer ns,i1,i2,nch,ib
1672       common/ipars/nch
1673       double precision sca(5,5,0:40000,2),sca1(5,5,0:40000,2)
1674       common/spac/sca,sca1
1675
1676       ns=900
1677       ksqma=8.2d0
1678       inck=ksqma/dble(ns)
1679       ksqmin=0.001d0
1680       lginck=(dlog(ksqma/ksqmin))/dble(ns)
1681
1682       print*,'Calculating screening amplitude'
1683
1684       do ib=0,ns+1
1685
1686          ksq=dble(ib-1)*inck
1687          lgksq=ksq
1688
1689          if(ib.eq.0)then
1690             ksq=0d0
1691             lgksq=0d0
1692          else
1693             lgksq=dble(ib-1)*lginck+dlog(ksqmin)
1694             ksq=dexp(lgksq)
1695          endif
1696
1697          do i1=1,nch
1698             do i2=1,nch
1699
1700                call screening(i1,i2,ksq,sc,sc1)
1701
1702                sca(i1,i2,ib,1)=lgksq
1703                sca(i1,i2,ib,2)=sc
1704                sca1(i1,i2,ib,1)=lgksq
1705                sca1(i1,i2,ib,2)=sc1
1706
1707             enddo
1708          enddo
1709
1710       enddo
1711
1712       return
1713       end
1714
1715 ccccc Interpolation....
1716
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)
1720       integer i,j,it
1721       common/opac/op,oph
1722
1723       incbt=op(1,1,2,1)-op(1,1,1,1)
1724       it=nint(bt/incbt)
1725       if(dble(it).gt.(bt/incbt))then
1726          it=it-1
1727       endif
1728
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)
1735
1736       out=m*del+op(i,j,it+1,2)
1737       out1=mh*delh+oph(i,j,it+1,2)
1738
1739       return
1740       end
1741
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)
1745       integer i,j,it,ns
1746       common/spac/sca,sca1
1747
1748       if(ktsq.lt.0.001d0)then
1749
1750          it=0
1751
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)
1758
1759          out=m*del+sca(i,j,it,2)
1760          out1=m1*del1+sca1(i,j,it,2)
1761
1762       elseif(ktsq.lt.8d0)then
1763
1764          ktmin=sca(1,1,1,1)
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
1768             it=it-1
1769          endif
1770
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)
1777
1778          out=m*del+sca(i,j,it+1,2)
1779          out1=m1*del1+sca1(i,j,it+1,2)
1780
1781       else
1782
1783          out=0d0
1784          out1=0d0
1785
1786       endif
1787
1788       return
1789       end
1790
1791 cccc  Integrates round Pomeron loop (to calculate screened amplitude)
1792
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
1801       common/ipars/nch
1802       common/spac/sca,sca1
1803       common/mom/q
1804       common/wvars/sig0,bb,t1,t2
1805       common/alphas/alphap,alpha0,alphapr,alpha0r,alphapm,alpha0m
1806
1807       do k=1,4
1808          a1(k)=q(k,1)-q(k,3)
1809          a2(k)=q(k,2)-q(k,4)
1810       enddo
1811
1812       nx=10
1813       hx=2d0/dble(nx)
1814
1815       out=0d0
1816
1817       do jx=-nx,nx
1818          do jy=-nx,nx
1819
1820             tpx=dble(jx)*hx
1821             tpy=dble(jy)*hx
1822             tp2=tpx**2+tpy**2
1823             wt=hx*hx
1824
1825             p1xp=p1x-tpx
1826             p1yp=p1y-tpy
1827             t12=p1xp**2+p1yp**2
1828             p2xp=tpx+p2x
1829             p2yp=tpy+p2y
1830             t22=p2xp**2+p2yp**2
1831
1832             a1(1)=-p1xp
1833             a1(2)=-p1yp
1834             a2(1)=-p2xp
1835             a2(2)=-p2yp
1836
1837       do i1=1,nch
1838          do i2=1,nch
1839
1840            call screeningint(i1,i2,tp2,sc,sc1)
1841            call wev(x0,a1,a2)
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))
1847             out=out+x0*wt*sc
1848
1849          enddo
1850       enddo
1851
1852
1853       enddo
1854       enddo
1855
1856       a1(1)=-p1x
1857       a1(2)=-p1y
1858       a2(1)=-p2x
1859       a2(2)=-p2y
1860
1861       t11=p1x**2+p1y**2
1862       t22=p2x**2+p2y**2
1863
1864       call formfac(t11,t22,x00p)
1865
1866       call wev(x00,a1,a2)
1867
1868
1869
1870       out=out+x00*x00p
1871
1872       return
1873       end
1874
1875 ccccc Pomeron -- diffrative eignstate i form factor
1876
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
1881       integer nch,i1,i2
1882       common/ipars/nch
1883
1884       out=0d0
1885
1886       delta1=dexp(-t1/2d0)
1887       delta2=dexp(-t2/2d0)
1888
1889       do i1=1,nch
1890          do i2=1,nch
1891
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))
1897
1898             out=out+wt
1899
1900          enddo
1901       enddo
1902
1903       return
1904       end
1905
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
1911       common/ipars/nch
1912       common/vars/s,rts,mmes,yx, iin
1913
1914       pi=dacos(-1d0)
1915
1916       nb=5000
1917       hb=99d0/dble(nb)
1918
1919       out=0d0
1920       out1=0d0
1921
1922       do ib=1,nb
1923          bt=dble(ib)*hb
1924          wt=-bt/2d0/pi*hb
1925
1926          call opacityint(i,j,bt,fr,fr1)
1927
1928          sige=sigo*dexp(dlog(rts)*2d0*ep)
1929          fr=fr*gaa(i)*gaa(j)*sige
1930          fr1=fr1*gaa(i)*gaa(j)*sige
1931
1932          out=out+wt*(1d0-dexp(-fr/2d0))*besj0(bt*dsqrt(ktsq))
1933      &        *gaa(i)*gaa(j)
1934          out1=out1+wt*(1d0-dexp(-fr1/2d0))*besj0(bt*dsqrt(ktsq))
1935      &        *gaa(i)*gaa(j)
1936
1937       enddo
1938
1939       return
1940       end
1941
1942
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
1948       common/ipars/nch
1949       common/vars/s,rts,mmes,yx, iin
1950
1951       pi=dacos(-1d0)
1952
1953       ampi=0.02d0
1954       amro=1d0
1955       a4=4d0*ampi
1956       alo=dlog(amro/ampi)
1957       bpol=2.4d0
1958
1959       nt=6000
1960       htt=6d0/dble(nt)**2
1961
1962          out=0d0
1963          out1=0d0
1964
1965          do it=0,nt
1966             t=dble(it)**2*htt
1967             if(it.eq.0) t=1d-8
1968             wt=htt*2d0*dble(it)/4d0/pi
1969             if(it.eq.0) wt=wt/2d0
1970             bes0=besj0(bt*dsqrt(t))
1971
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))
1976
1977             asp1=asp
1978             form1=dlog(ffi*ffj)-2d0*t*asp1*dlog(rts)
1979 cccccccccccccccc
1980             ah=dsqrt(1d0+a4/t)
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)
1983 ccccccccccccccc
1984             ww=bes0*dexp(form1-2d0*h1pi*dlog(rts))
1985             aspt=t*asp+h1pi
1986
1987             out=out+ww*wt
1988             out1=out1+bes0*wt*ffi*ffj
1989
1990          enddo
1991
1992
1993       return
1994       end
1995
1996       function rgauss(r1,r2)
1997       implicit double precision (a-y)
1998
1999       pi=dacos(-1d0)
2000       rgauss=dsqrt(-2d0*dlog(r1))*dsin(2d0*pi*r2)
2001
2002       return
2003       end
2004
2005 c     boosting subroutine
2006
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))
2012          do 10 j=1,3
2013   10     plb(j)=pcm(j)+fact*pboo(j)
2014       return
2015       end
2016
2017 c     calculates lorentzian dot product for real 4-vectors
2018
2019       subroutine dot(v1,v2,dt)
2020       double precision v1(4),v2(4),dt
2021
2022       dt=-(v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3)
2023      &-v1(4)*v2(4))
2024
2025       return
2026       end
2027 c
2028
2029 *
2030 * subtractive mitchell-moore generator
2031 * ronald kleiss - october 2, 1987
2032 *
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)
2042 *
2043 * in this version  M =2**30  and  RM=1/M=2.D0**(-30.D0)
2044 *
2045 * the array NR has been initialized by putting NR(i)=i
2046 * and subsequently running the algorithm 100,000 times.
2047 *
2048
2049       subroutine R2455(Ran)
2050       IMPLICIT REAL*8(A-H,O-Z)
2051       DIMENSION N(55)
2052       DATA N/
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/
2064       DATA M/1073741824/
2065       DATA RM/0.9313225746154785D-09/
2066       DATA K/55/,L/31/
2067       IF(K.EQ.55) THEN
2068          K=1
2069       ELSE
2070          K=K+1
2071       ENDIF
2072       IF(L.EQ.55) THEN
2073          L=1
2074       ELSE
2075          L=L+1
2076       ENDIF
2077       J=N(L)-N(K)
2078       IF(J.LT.0) J=J+M
2079       N(K)=J
2080       RAN=J*RM
2081       END
2082
2083
2084
2085       SUBROUTINE RAMBO(N,ET,XM,P,WT)
2086 C------------------------------------------------------
2087 C
2088 C                       RAMBO
2089 C
2090 C    RA(NDOM)  M(OMENTA)  B(EAUTIFULLY)  O(RGANIZED)
2091 C
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
2095 C
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
2101 C
2102 C------------------------------------------------------
2103       implicit none
2104 !     boost arguments
2105       integer n
2106       double precision et,xm(100),p(4,100),wt
2107 !     function
2108       double precision rn
2109 !     internal variables
2110       double precision q(4,100),z(100),r(4),
2111      &     b(3),p2(100),xm2(100),e(100),v(100)
2112       integer iwarn(5)
2113       double precision acc
2114       integer itmax,ibegin
2115       data acc/1.0d-14/,itmax/6/,ibegin/0/,iwarn/5*0/
2116
2117       integer i,k,iter
2118       double precision twopi,po2log,xmt,c,s,f,rmas,g,a,x,bq,
2119      &     accu,f0,g0,wt2,wt3,wtm,x2,xmax
2120       integer nm
2121       save
2122
2123
2124 C
2125 C INITIALIZATION STEP: FACTORIALS FOR THE PHASE SPACE WEIGHT
2126       IF(IBEGIN.NE.0) GOTO 103
2127       IBEGIN=1
2128       TWOPI=8.*DATAN(1.D0)
2129       PO2LOG=DLOG(TWOPI/4.)
2130       Z(2)=PO2LOG
2131       DO 101 K=3,100
2132   101 Z(K)=Z(K-1)+PO2LOG-2.*DLOG(DFLOAT(K-2))
2133       DO 102 K=3,100
2134   102 Z(K)=(Z(K)-DLOG(DFLOAT(K-1)))
2135 C
2136 C CHECK ON THE NUMBER OF PARTICLES
2137   103 IF(N.GT.1.AND.N.LT.101) GOTO 104
2138       PRINT 1001,N
2139       STOP
2140 C
2141 C CHECK WHETHER TOTAL ENERGY IS SUFFICIENT; COUNT NONZERO MASSES
2142   104 XMT=0.
2143       NM=0
2144       DO 105 I=1,N
2145       IF(XM(I).NE.0.D0) NM=NM+1
2146   105 XMT=XMT+DABS(XM(I))
2147       IF(XMT.LE.ET) GOTO 201
2148       PRINT 1002,XMT,ET
2149       STOP
2150 C
2151 C THE PARAMETER VALUES ARE NOW ACCEPTED
2152 C
2153 C GENERATE N MASSLESS MOMENTA IN INFINITE PHASE SPACE
2154   201 DO 202 I=1,N
2155       C=2.*RN(1)-1.
2156       S=DSQRT(1.-C*C)
2157       F=TWOPI*RN(2)
2158       Q(4,I)=-DLOG(RN(3)*RN(4))
2159       Q(3,I)=Q(4,I)*C
2160       Q(2,I)=Q(4,I)*S*DCOS(F)
2161  202  Q(1,I)=Q(4,I)*S*DSIN(F)
2162
2163 C
2164 C CALCULATE THE PARAMETERS OF THE CONFORMAL TRANSFORMATION
2165       DO 203 I=1,4
2166   203 R(I)=0.
2167       DO 204 I=1,N
2168       DO 204 K=1,4
2169   204 R(K)=R(K)+Q(K,I)
2170       RMAS=DSQRT(R(4)**2-R(3)**2-R(2)**2-R(1)**2)
2171       DO 205 K=1,3
2172   205 B(K)=-R(K)/RMAS
2173       G=R(4)/RMAS
2174       A=1./(1.+G)
2175       X=ET/RMAS
2176 C
2177 C TRANSFORM THE Q'S CONFORMALLY INTO THE P'S
2178       DO 207 I=1,N
2179       BQ=B(1)*Q(1,I)+B(2)*Q(2,I)+B(3)*Q(3,I)
2180       DO 206 K=1,3
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)
2183 C
2184 C CALCULATE WEIGHT AND POSSIBLE WARNINGS
2185       WT=PO2LOG
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
2189       IWARN(1)=IWARN(1)+1
2190   208 IF(WT.LE. 174.D0) GOTO 209
2191       IF(IWARN(2).LE.5) PRINT 1005,WT
2192       IWARN(2)=IWARN(2)+1
2193 C
2194 C RETURN FOR WEIGHTED MASSLESS MOMENTA
2195   209 IF(NM.NE.0) GOTO 210
2196       WT=DEXP(WT)
2197       RETURN
2198 C
2199 C MASSIVE PARTICLES: RESCALE THE MOMENTA BY A FACTOR X
2200   210 XMAX=DSQRT(1.-(XMT/ET)**2)
2201       DO 301 I=1,N
2202       XM2(I)=XM(I)**2
2203   301 P2(I)=P(4,I)**2
2204       ITER=0
2205       X=XMAX
2206       ACCU=ET*ACC
2207   302 F0=-ET
2208       G0=0.
2209       X2=X*X
2210       DO 303 I=1,N
2211       E(I)=DSQRT(XM2(I)+X2*P2(I))
2212       F0=F0+E(I)
2213   303 G0=G0+P2(I)/E(I)
2214       IF(DABS(F0).LE.ACCU) GOTO 305
2215       ITER=ITER+1
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)
2221       ENDIF
2222       GOTO 305
2223   304 X=X-F0/(X*G0)
2224       GOTO 302
2225   305 DO 307 I=1,N
2226       V(I)=X*P(4,I)
2227       DO 306 K=1,3
2228   306 P(K,I)=X*P(K,I)
2229   307 P(4,I)=E(I)
2230 C
2231 C CALCULATE THE MASS-EFFECT WEIGHT FACTOR
2232       WT2=1.
2233       WT3=0.
2234       DO 308 I=1,N
2235       WT2=WT2*V(I)/E(I)
2236   308 WT3=WT3+V(I)**2/E(I)
2237       WTM=(2.*N-3.)*DLOG(X)+DLOG(WT2/WT3*ET)
2238 C
2239 C RETURN FOR  WEIGHTED MASSIVE MOMENTA
2240       WT=WT+WTM
2241       IF(WT.GE.-180.D0) GOTO 309
2242       IF(IWARN(3).LE.5) PRINT 1004,WT
2243       IWARN(3)=IWARN(3)+1
2244   309 IF(WT.LE. 174.D0) GOTO 310
2245       IF(IWARN(4).LE.5) PRINT 1005,WT
2246       IWARN(4)=IWARN(4)+1
2247   310 WT=DEXP(WT)
2248       RETURN
2249 C
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)
2258       END
2259 C
2260
2261       function rn(idum)
2262 *
2263 * SUBTRACTIVE MITCHELL-MOORE GENERATOR
2264 * RONALD KLEISS - OCTOBER 2, 1987
2265 *
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)
2275 *
2276 * IN THIS VERSION  M =2**30  AND  RM=1/M=2.D0**(-30D0)
2277 *
2278 * THE ARRAY NR HAS BEEN INITIALIZED BY PUTTING NR(I)=I
2279 * AND SUBSEQUENTLY RUNNING THE ALGORITHM 100,000 TIMES.
2280 *
2281
2282       implicit none
2283       double precision rn
2284       integer idum
2285       integer n(55)
2286       data n/
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
2299       double precision rm
2300       integer j,k,l,m
2301
2302       data eps/1D-9/
2303       data m/1073741824/
2304       data rm/0.9313225746154785D-09/
2305       data k/55/,l/31/
2306
2307
2308   1   CONTINUE
2309       IF(K.EQ.55) THEN
2310          K=1
2311       ELSE
2312          K=K+1
2313       ENDIF
2314       IF(L.EQ.55) THEN
2315          L=1
2316       ELSE
2317          L=L+1
2318       ENDIF
2319       J=N(L)-N(K)
2320       IF(J.LT.0) J=J+M
2321       N(K)=J
2322       RN=J*RM
2323       IF(RN.LT.EPS) GOTO 1
2324       IF(RN.GT.1D0-EPS) GOTO 1
2325       RETURN
2326       END