EVGEN/PartEloss/integralgauss.f

   1 *******************************************************************
   2       SUBROUTINE eloss_qsimp(sr,sd)
   3       REAL*8           xr,length,noc,qq,nnorm,rri,fracc
   4       COMMON /input/   xr,length,noc,qq,nnorm,rri,fracc
   5       INTEGER JMAX
   6       REAL*8 a1,b1,diff1,d1,sr,sd,EPS,fu1r,fu1d
   7       EXTERNAL eloss_func1, eloss_func2
   8       PARAMETER (EPS=1.e-9, JMAX=12)
   9       REAL*8 osr,ostr
  10       ostr=-1.e10
  11       osr=-1.e10
  12       a1=0.0
  13       b1=(0.99999-xr)*fracc
  14       diff1 = b1-a1
  15       d1 = 0.5*diff1
  16 *
  17       res=dgauss(eloss_func2,a1,b1,1.d-6)
  18       call eloss_func1(a1,fu1r,fu1d)
  19       sd=fu1d
  20       sr=res
  21       END
  22
  23       function eloss_func2(yy)
  24       implicit double precision (a-h,o-z)
  25       call eloss_func1(yy,fu1r,fu1d)
  26       eloss_func2=fu1r
  27       return
  28       end
  29
  30
  31 **************************************************************
  32 *
  33       SUBROUTINE eloss_func1(yy,funr,fund)
  34 *
  35       REAL*8           funr,yy,fund
  36       REAL*8           xr,length,noc,qq,nnorm,rri,fracc
  37       COMMON /input/   xr,length,noc,qq,nnorm,rri,fracc
  38       EXTERNAL         eloss_lookup
  39       REAL*8           cont, disc, wwt, tepsi
  40 *
  41       tepsi = yy
  42       wwt = tepsi
  43       if(wwt.ge.1.3) then
  44          call eloss_lookup(rri,1.d0,cont,disc)
  45          funr=0.0
  46          fund=disc
  47       else
  48          call eloss_lookup(rri,wwt,cont,disc)
  49          funr = cont*eloss_fragm(xr/(1.0-tepsi/fracc),qq)
  50      .         /(1.0-tepsi/fracc)
  51          fund = disc
  52       endif
  53       END
  54 *******************************************************************
  55       SUBROUTINE eloss_lookup(rrrr,xxxx,continuous,discrete)
  56 *
  57       REAL*8           xx(400), da(30), ca(30,260), rrr(30)
  58       COMMON /data/    xx, da, ca, rrr
  59       REAL*8           rrrr,xxxx, continuous, discrete
  60       REAL*8           rrin, xxin
  61       INTEGER          nrlow, nrhigh, nxlow, nxhigh
  62       REAL*8           rrhigh, rrlow, rfraclow, rfrachigh
  63       REAL*8           xfraclow, xfrachigh
  64       REAL*8           clow, chigh
  65 *
  66       rrin = rrrr
  67       xxin = xxxx
  68 *
  69 *    determine the tabulated values xx(nxlow), xx(nxhigh)
  70 *    rrlow, rrhigh such that
  71 *    xx(nxlow) < xxin <  xx(nxhigh)
  72 *    rrlow < rrin < rrhigh
  73 *
  74       nxlow = int(xxin/0.005) + 1
  75       nxhigh = nxlow + 1
  76       xfraclow = (xx(nxhigh)-xxin)/0.005
  77       xfrachigh = (xxin - xx(nxlow))/0.005
  78 *
  79       do 666, nr=1,30
  80          if (rrin.lt.rrr(nr)) then
  81             rrhigh = rrr(nr)
  82          else
  83             rrhigh = rrr(nr-1)
  84             rrlow = rrr(nr)
  85             nrlow = nr
  86             nrhigh = nr-1
  87             goto 665
  88          endif
  89  666     enddo
  90  665     continue
  91 *
  92       rfraclow = (rrhigh-rrin)/(rrhigh-rrlow)
  93       rfrachigh = (rrin-rrlow)/(rrhigh-rrlow)
  94 *
  95       clow = xfraclow*ca(nrlow,nxlow)+xfrachigh*ca(nrlow,nxhigh)
  96       chigh = xfraclow*ca(nrhigh,nxlow)+xfrachigh*ca(nrhigh,nxhigh)
  97       continuous = rfraclow*clow + rfrachigh*chigh
  98       discrete = rfraclow*da(nrlow) + rfrachigh*da(nrhigh)
  99 *
 100       END
 101
 102 ***************************************************************
 103
 104 c       BKK FF
 105
 106       FUNCTION eloss_fragmbkk(xxx,qqq)
 107       REAL*8   alphav, betav, gammav, nv
 108       REAL*8   alphas, betas, gammas, ns
 109       REAL*8   sbar, xx, qq, xxx, qqq, lambda, fragv, frags
 110 *
 111       xx = xxx
 112       qq = qqq
 113       lambda = 1.0
 114       sbar=log(log(qq*qq/(lambda*lambda))/log(4.0/(lambda*lambda)))
 115 *
 116       alphav = -1.0 - 0.0272*sbar
 117       betav = 1.2 + 0.67*sbar
 118       gammav = -0.393*sbar
 119       nv = 0.551 - 0.053*sbar - 0.032*sbar*sbar
 120 *
 121 *      alphav = -1.0 - 0.059*sbar
 122 *      betav = 1.2 + 0.6*sbar
 123 *      gammav = -0.163*sbar
 124 *      nv = 0.338 - 0.064*sbar - 0.0105*sbar*sbar
 125 *
 126       alphas = -1.0 + 0.447*sbar - 0.266*sbar*sbar
 127       betas = 4.7 - 2.88*sbar + 2.05*sbar*sbar
 128       gammas = -9.01*sbar + 4.36*sbar*sbar
 129       ns = 1.23 + 2.85*sbar - 1.6*sbar*sbar
 130 *
 131 *      alphas = -1.0 + 0.757*sbar - 0.537*sbar*sbar
 132 *      betas = 5.26 - 5.22*sbar + 3.62*sbar*sbar
 133 *      gammas = -13.6*sbar + 8.17*sbar*sbar
 134 *      ns = 1.19 + 4.20*sbar - 2.86*sbar*sbar
 135 *
 136       fragv = nv*(xx**alphav)*((1.0-xx)**betav)*((1.0+xx)**gammav)
 137       frags = ns*(xx**alphas)*((1.0-xx)**betas)*((1.0+xx)**gammas)
 138 c      fragmbkk = (fragv+frags)
 139       eloss_fragmbkk = fragv
 140       END
 141
 142 ***************************************************************
 143
 144       FUNCTION eloss_fragm(xxx,qqq)
 145       REAL*8   alpha, beta, gamma, n
 146       REAL*8   sbar, xx, qq, xxx, qqq, lambda
 147 *
 148       xx = xxx
 149       qq = qqq
 150       lambda = 0.088
 151       sbar=log(log(qq*qq/(lambda*lambda))/log(2.0/(lambda*lambda)))
 152 *
 153
 154 c       u or d -> pi+ + pi-
 155
 156       n=0.54610-0.22946*sbar-0.22594*sbar**2+0.21119*sbar**3
 157       alpha=-1.46616-0.45404*sbar-0.12684*sbar**2+0.27646*sbar**3
 158       beta=1.01864+0.95367*sbar-1.09835*sbar**2+0.74657*sbar**3
 159       gamma=-0.01877*sbar+0.02949*sbar**2
 160
 161 c       g -> pi+ + pi-
 162
 163       n=6.04510-6.61523*sbar-1.64978*sbar**2+2.68223*sbar**3
 164       alpha=-.71378+0.14705*sbar-1.08423*sbar**2-.43182*sbar**3
 165       beta=2.92133+1.48429*sbar+1.32887*sbar**2-1.78696*sbar**3
 166       gamma=0.23086*sbar-0.29182*sbar**2
 167
 168       eloss_fragm = n*xx**alpha*(1.-xx)**beta*(1.+gamma/xx)/2.
 169
 170       END
 171