TUHKMgen/PYQUEN/progs_fortran.f

   1 c*******************************************
   2
   3       SUBROUTINE myini
   4
   5       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
   6       double precision nbcol,npart,npart0
   7       real p,v,plu
   8       external ludata,pydata
   9       double precision numpar,npar0,nbco0
  10 c      common /lujets/ n,k(150000,5),p(150000,5),v(150000,5)
  11       common /hyjets/ nhj,nhp,khj(150000,5),phj(150000,5),vhj(150000,5)
  12       common /hyfpar/ bgen,nbcol,npart,npart0,npyt,nhyd
  13       common /hyflow/ ytfl,ylfl,Tf,fpart
  14       common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet
  15       common /pydat1/ mstu(200),paru(200),mstj(200),parj(200)
  16       common /pysubs/ msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
  17       common /pypars/ mstp(200),parp(200),msti(200),pari(200)
  18       common /pyqpar/ T0,tau0,nf,ienglu,ianglu
  19       common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
  20       common /hypyin/ ene,rzta,rnta,bfix,ifb,nh
  21
  22       COMMON/PYDATR/MRPY(6),RRPY(100)
  23       common/SERVICE/iseed_fromC,iPythDecay,parPYTH(100)
  24
  25
  26       common /hypart/ppart(10,50000),bmin,bmax,njp
  27 c      save /lujets/,/hyjets/,/hyflow/,/hyjpar/,/hyfpar/,/pyqpar/,
  28 c     >     /pysubs/,/pypars/,/pydat1/
  29
  30 * ----------- INITIALIZATION OF RUNDOM GENERATOR
  31       MRPY(1)=iseed_fromC
  32       write(*,*)'--initialization of high-pt part-- :  '
  33 *      write(*,*)'evnt generator phase  ',iseed_fromC
  34
  35 * initialize HYINIT with the input parameters
  36       write(*,*)' AW= ', AW,' energy= ', ene, ' ptmin= ',ptmin
  37       write(*,*)' bminh= ',bminh,' bmaxh= ', bmaxh
  38       write(*,*)' ifb= ', ifb, ' bfix= ', bfix
  39       write(*,*)' ishad= ',ishad, ' nhsel=  ',nhsel
  40       write(*,*)' ienglu= ',ienglu, ' ianglu= ',ianglu
  41       write(*,*)' T0= ',T0,' tau0= ',tau0,' nf= ', nf
  42
  43 * set input PYTHIA parameters from SERVICE common:
  44       CKIN(3)=ptmin               ! minimum pt in initial hard scattering, GeV
  45       MSEL=0                      ! QCD-dijet production
  46       MSTU(21)=1                  ! avoid stopping run
  47       PARU(14)=1.              ! tolerance parameter to adjust fragmentation
  48       MSTP(81)=1                  ! pp multiple scattering off
  49       MSTJ(21)=0                    ! hadron decays off
  50       MSTP(2)=1                !which order running alphaS
  51       MSTP(33)=0               !unclusion of k factor in cross section
  52       MSTP(51)=7               !structure function chosen - CTEQ5M pdf
  53       MSTP(82)=4               !Defines the multi-parton model
  54       PARP(67)=1              !amount of initial-state radiation
  55       PARP(82)=2.              !pt cutoff for multiparton interactions
  56       MSTJ(11)=3               !Choice of the fragmentation function
  57
  58 * flags to select dijet production channels (active if msel=0)
  59       MSUB(11)=1   ! q+q->q+q
  60       MSUB(12)=1   ! q+qbar->q+qbar
  61       MSUB(13)=1   ! q+qbar->g+g
  62       MSUB(28)=1   ! q+g->q+g
  63       MSUB(53)=1   ! g+g->q+qbar
  64       MSUB(68)=1   ! g+g->g+g
  65
  66 * flags to select prompt photon production channels (active if msel=0)
  67       MSUB(14)=1   ! q+qbar->g+gamma
  68       MSUB(18)=1   ! q+qbar->gamma+gamma
  69       MSUB(29)=1   ! q+g->q+gamma
  70       MSUB(114)=1  ! g+g->gamma+gamma
  71       MSUB(115)=1  ! g+g->g+gamma
  72
  73
  74       write(*,*)'--initialization of PYTHIA parameters-- :  '
  75
  76       write(*,*)'CKIN(3) MSEL MSTU(21) PARU(14)'
  77       write(*,*)CKIN(3), MSEL, MSTU(21), PARU(14)
  78
  79       write(*,*)'MSTP(81) MSTJ(21) MSTP(2) MSTP(33)'
  80       write(*,*)MSTP(81), MSTJ(21), MSTP(2), MSTP(33)
  81
  82       write(*,*)'MSTP(51) MSTP(82) PARP(67) PARP(82) MSTJ(11)'
  83       write(*,*)MSTP(51),MSTP(82),PARP(67),PARP(82),MSTJ(11)
  84
  85       write(*,*)'MSUB(11) MSUB(12) MSUB(13) MSUB(28) MSUB(53) MSUB(68)'
  86       write(*,*)MSUB(11),MSUB(12),MSUB(13),MSUB(28),MSUB(53),MSUB(68)
  87
  88       write(*,*)'MSUB(14) MSUB(18) MSUB(29) MSUB(114) MSUB(115)'
  89       write(*,*)MSUB(14),MSUB(18),MSUB(29),MSUB(114),MSUB(115)
  90
  91       call hyinit     !       (ene,AW,ifb,bmin,bmax,bfix1,nh)
  92
  93       end
  94 ****************end myini***************************************************************
  95
  96
  97
  98 ********************************* HYINIT ***************************
  99       SUBROUTINE HYINIT
 100
 101 *ml--      (energy,A,ifb1,bmin,bmax,bfix1,nh1)
 102 *     PYTHIA inizialization, calculation of total and hard cross sections and
 103 *     # of participants and binary sub-collisions at "reference point" (Pb,b=0),
 104 *     tabulation of nuclear thickness function and nuclear overlap function,
 105 *     test of input parametes
 106       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
 107       double precision numpar,npar0,nbco0
 108       external numpar,rhoaa,hfunc3
 109       common /pyint7/ sigt(0:6,0:6,0:5)
 110       common /pypars/ mstp(200),parp(200),msti(200),pari(200)
 111       common /pysubs/ msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
 112       common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet
 113       common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
 114       common /hypyin/ ene,rzta,rnta,bfix,ifb,nh
 115       common /hyflow/ ytfl,ylfl,Tf,fpart
 116       common /hygeom/ BC
 117       common /hythic/ BAB(110),TAB(110),TAAB(110)
 118       common /hynup1/ bp,x
 119
 120 *
 121
 122 c--      save /pyint7/,/pypars/,/pysubs/,/hyjpar/,/hyipar/,/hypyin/,
 123 c--     >     /hyflow/,/hygeom/,/hythic/
 124
 125 * start HYDJET initialization
 126       init=1
 127       ipr=1
 128
 129 * set beam paramters
 130 c-ml      ene=energy                               ! c.m.s. energy per nucleon
 131 c-ml      AW=A                                     ! atomic weight
 132
 133
 134       RA=1.15d0*AW**0.333333d0                 ! nuclear radius
 135       rzta=1.d0/(1.98d0+0.015d0*A**0.666667d0) ! fraction of protons in nucleus
 136       rnta=1.d0-rzta                           ! fraction of neutrons in nucleus
 137 c-ml      ifb=ifb1                                 ! centrality flag
 138 c-ml      bfix=bfix1                               ! fixed impact parameter
 139 c-ml      nh=nh1                                   ! mean soft mult. in central PbPb
 140
 141 c-ml      ptmin=ckin(3)                            ! minimum pt of hard scattering
 142
 143 c* printout of HYDJET versioning
 144 c      write(6,*) '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 145 c     >%%%%%%%%%%%%%%%%%%%%%%%%'
 146 c      write(6,*) '%%                        This is HYDJET version 1.5
 147 c     >                      %%'
 148 c      write(6,*) '%%         Corresponding author: Igor Lokhtin (Igor.Lo
 149 c     >khtin@cern.ch)        %%'
 150 c      write(6,*) '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 151 c     >%%%%%%%%%%%%%%%%%%%%%%%%'
 152
 153 * Pythia inizialization for pp collisions
 154       call pyinit('cms','p','p',ene)
 155
 156 * calculation of hard scattering cross section
 157       if(nhsel.ne.0) then
 158
 159        mstp(111)=0
 160 * no printout of Pythia initialization information hereinafter
 161        mstp(122)=0
 162
 163
 164 * initialize HYINIT with the input parameters
 165       write(*,*)' in hyinit AW= ', AW,' energy= ', ene, ' ptmin= ',ptmin
 166       write(*,*)' bminh= ',bminh,' bmaxh= ', bmaxh
 167       write(*,*)' ifb= ', ifb, ' bfix= ', bfix
 168
 169
 170 * Pythia test pp event run
 171        do i=1,1000
 172         call pyevnt
 173        end do
 174 * hard scattering pp cross section
 175        sjpp=pari(1)
 176
 177 * Pythia inizialization for pn collisions
 178       call pyinit('cms','p','n',ene)
 179 * Pythia test pn event run
 180        do i=1,1000
 181         call pyevnt
 182        end do
 183 * hard scattering pn cross section
 184        sjpn=pari(1)
 185
 186 * Pythia inizialization for nn collisions
 187       call pyinit('cms','n','n',ene)
 188 * Pythia test nn event run
 189        do i=1,1000
 190         call pyevnt
 191        end do
 192 * hard scattering nn cross section
 193        sjnn=pari(1)
 194
 195        sigjet=rzta*rzta*sjpp+rnta*rnta*sjnn+2.d0*rzta*rnta*sjpn
 196
 197       end if
 198
 199 * total inelastic cross section
 200       if(sigin.lt.10.d0.or.sigin.gt.200.d0)
 201      >   sigin=sigt(0,0,0)-sigt(0,0,1)
 202
 203 * # of nucelons-participants and NN sub-collisions at "reference point" (Pb,b=0)
 204       Apb=207.d0
 205       Rpb=1.15d0*Apb**0.333333d0
 206       EPS=0.005d0
 207       Z2=4.d0*Rpb
 208       Z1=-1.d0*Z2
 209       H=0.01d0*(Z2-Z1)
 210       do ib=1,110
 211        BC=3.d0*Rpb*(ib-1)/109.d0
 212        CALL SIMPA(Z1,Z2,H,EPS,1.d-8,rhoaa,Z,RES,AIH,AIABS)
 213        BAB(ib)=BC
 214        TAB(ib)=Apb*RES
 215       end do
 216       Z1=0.d0
 217       Z2=6.28318d0
 218       H=0.01d0*(Z2-Z1)
 219       bp=0.d0
 220       CALL SIMPA(Z1,Z2,H,EPS,1.d-8,HFUNC3,X,TAAPB0,AIH,AIABS)
 221
 222       npar0=numpar(0.d0)                            ! Npart(Pb,b=0)
 223       nbco0=0.1d0*sigin*TAAPB0                      ! Nsub(Pb,b=0)
 224
 225       init=0
 226
 227 * creation of arrays for tabulation of beam/target nuclear thickness function
 228       Z2=4.d0*RA
 229       Z1=-1.d0*Z2
 230       H=0.01d0*(Z2-Z1)
 231       do ib=1,110
 232        BC=3.d0*RA*(ib-1)/109.d0
 233        CALL SIMPA(Z1,Z2,H,EPS,1.d-8,rhoaa,Z,RES,AIH,AIABS)
 234        BAB(ib)=BC
 235        TAB(ib)=AW*RES
 236       end do
 237
 238 * creation of arrays for tabulation of nuclear overlap function
 239       Z1=0.d0
 240       Z2=6.28318d0
 241       H=0.01d0*(Z2-Z1)
 242       do ib=1,110
 243        bp=BAB(ib)
 244        CALL SIMPA(Z1,Z2,H,EPS,1.d-8,HFUNC3,X,RES,AIH,AIABS)
 245        TAAB(ib)=RES
 246       end do
 247
 248 * test of centrality selection
 249       if(ifb.eq.0) then
 250        if(bfix.lt.0.d0) then
 251         write(6,*) 'Impact parameter less than zero!'
 252         bfix=0.d0
 253        end if
 254        if (bfix.gt.3.d0) then
 255         write(6,*) 'Impact parameter larger than three nuclear radius!'
 256         bfix=3.d0
 257        end if
 258       else
 259        if(bmin.lt.0.d0) then
 260         write(6,*) 'Impact parameter less than zero!'
 261         bmin=0.d0
 262        end if
 263        if(bmax.gt.3.d0) then
 264         write(6,*) 'Impact parameter larger than three nuclear radius!'
 265         bmax=3.d0
 266        end if
 267 c-ml       bminh=bmin
 268 c-ml       bmaxh=bmax
 269       end if
 270
 271 * test of flow parameter selection
 272       if (Tf.lt.0.08d0.or.Tf.gt.0.2d0) Tf=0.1d0       ! freeze-out temperature
 273       if (ylfl.lt.0.01d0.or.ylfl.gt.7.d0) ylfl=4.d0   ! longitudinal flow rapidity
 274       if (ytfl.lt.0.01d0.or.ytfl.gt.3.d0) ytfl=1.5d0  ! transverse flow rapidity
 275       if (fpart.le.0.d0.or.fpart.gt.1.d0) fpart=1.d0  ! fraction of soft multiplicity
 276                                                       ! proport. to # of participants
 277 * test of 'nhsel' selection
 278       if(nhsel.ne.1.and.nhsel.ne.2.and.nhsel.ne.3.and.nhsel.ne.4)
 279      > nhsel=0
 280
 281       return
 282       end
 283 ****************************** END HYINIT ***************************
 284       SUBROUTINE MYDELTA
 285       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
 286       common/SERVICEEV/delta,KC,ipdg
 287       COMMON /PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
 288       integer KC, PYCOMP
 289       real delta
 290 * allows to transmit delta -limits for BW from PYTHIA to c++ part
 291        KC=PYCOMP(ipdg)
 292        delta=PMAS(int(KC),3)
 293 c       write(*,*)" ipdg ", ipdg, " KC ",KC," delta ",delta
 294        return
 295        end
 296
 297
 298 ********************************* HYEVNT ****************************
 299       SUBROUTINE HYEVNT
 300 *     generation of single HYDJET event (soft+hard parts) at given parameters
 301       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
 302       double precision numpar,npar0,nbco0,npart,nbcol,npart0
 303       real p,v, delta
 304       external hsin,hftaa,numpar,hyhard,hipsear,pyr,pymass,PYCOMP
 305       external ludata
 306 c      common /lujets/ n,k(150000,5),p(150000,5),v(150000,5)
 307       common /hyjets/ nhj,nhp,khj(150000,5),phj(150000,5),vhj(150000,5)
 308       common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
 309       common /hypyin/ ene,rzta,rnta,bfix,ifb,nh
 310       common /hyfpar/ bgen,nbcol,npart,npart0,npyt,nhyd
 311       common /hyflow/ ytfl,ylfl,Tf,fpart
 312       common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet
 313 c-ml
 314       common /hypart/ppart(10,50000),bmin,bmax,njp
 315       COMMON /PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
 316
 317 c--      save/lujets/,/hyjets/,/hyipar/,/hyfpar/,/hyflow/,/hyjpar/,/hypyin/
 318
 319
 320 * reset lujets and hyjets arrays before event generation
 321
 322          write(*,*)'in hyevnt 0'
 323
 324       n=0
 325       nhj=0
 326
 327       do ncl=1,150000
 328        do j=1,5
 329 c        p(ncl,j)=0.
 330         phj(ncl,j)=0.d0
 331 c        v(ncl,j)=0.
 332         vhj(ncl,j)=0.d0
 333 c        k(ncl,j)=0
 334         khj(ncl,j)=0
 335        enddo
 336       end do
 337 *
 338       pi=3.14159d0
 339
 340 * generate impact parameter of A-A collision
 341       if(ifb.eq.0) then
 342        b1=bfix*RA
 343        bgen=bfix
 344       else
 345        call hipsear(fmax1,xmin1)
 346        fmax=fmax1
 347        xmin=xmin1
 348  3     bb1=xmin*pyr(0)+bminh*RA
 349        ff1=fmax*pyr(0)
 350        fb=hsin(bb1)
 351        if(ff1.gt.fb) goto 3
 352        b1=bb1
 353        bgen=bb1/RA
 354       end if
 355
 356        write(*,*)'in hyevnt bminh bgen sigin',RA,bminh,bgen,sigin
 357
 358 * calculate # of nucelons-participants and binary NN sub-collisions
 359       npart=numpar(b1)                          ! Npart(b)
 360       npart0=numpar(0.d0)                          ! Npart(b)
 361
 362 c      write(*,*)'in hyevnt npart',npart
 363       nbcol=0.1d0*sigin*hftaa(b1)               ! Nsub(b)
 364 c      write(*,*)'in hyevnt nbcol',nbcol
 365 * generate hard parton-parton scatterings (Q>ptmin) 'njet' times
 366       njet=0
 367       if(nhsel.ne.0) then
 368        pjet=sigjet/sigin
 369 c       write(*,*)'in hyevnt pjet',pjet
 370        do i=1,int(nbcol)
 371         if(pyr(0).lt.pjet) njet=njet+1
 372        end do
 373 c       write(*,*)'before hyhard'
 374        call hyhard
 375       end if
 376
 377       npyt=nhj
 378
 379 c      write(*,*)'in hyevnt pjet njet=', pjet,njet
 380 c-ml
 381 * fill array 'ppart'
 382
 383       do ih=1,nhj
 384         if(i.le.50000)then
 385         ppart(1,ih)= khj(ih,1) ! pdg ?
 386         ppart(2,ih)= khj(ih,2) ! pdg  -mother
 387         ppart(3,ih)=phj(ih,1) ! px
 388         ppart(4,ih)=phj(ih,2) ! py
 389         ppart(5,ih)=phj(ih,3) ! pz
 390         ppart(6,ih)=phj(ih,4) ! E
 391         ppart(7,ih)= vhj(ih,1)
 392         ppart(8,ih)= vhj(ih,2)
 393         ppart(9,ih)= vhj(ih,3)
 394         ppart(10,ih)= vhj(ih,4)
 395
 396         endif
 397       end do
 398
 399        njp=nhj ! fill number of jet particles for InitialStateBjorken
 400 c--
 401
 402 * fill array 'lujets'
 403
 404 c      do ih=1,nhj
 405 c       do jh=1,5
 406 c        p(ih,jh)=phj(ih,jh)
 407 c        k(ih,jh)=khj(ih,jh)
 408 c        v(ih,jh)=vhj(ih,jh)
 409 c       end do
 410 c      end do
 411
 412 c-ml
 413 c-      call luedit(2)                      ! remove unstable particles and partons
 414
 415
 416
 417       return
 418       end
 419 ****************************** END HYEVNT **************************
 420 ********************************* HYHARD ***************************
 421       SUBROUTINE HYHARD
 422 *     generate 'njet' number of hard parton-parton scatterings with PYTHIA
 423       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
 424       double precision npar0,nbco0,npart,nbcol,npart0
 425       double precision gaus
 426
 427       INTEGER PYK,PYCHGE,PYCOMP
 428       CHARACTER beam*2,targ*2
 429       external pydata
 430       external pyp,pyr,pyk,pyquen,shad1,gaus
 431
 432       common /pyjets/ n,npad,k(4000,5),p(4000,5),v(4000,5)
 433       common /hyjets/ nhj,nhp,khj(150000,5),phj(150000,5),vhj(150000,5)
 434       COMMON /PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
 435       COMMON /PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
 436       COMMON /PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
 437       COMMON /PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
 438       common /pypars/ mstp(200),parp(200),msti(200),pari(200)
 439       common /parimp/ b1, psib1, r0, rb1, rb2, noquen
 440       common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet
 441       common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
 442       common /hyfpar/ bgen,nbcol,npart,npart0,npyt,nhyd
 443       common /hypyin/ ene,rzta,rnta,bfix,ifb,nh
 444 c--      save /pyjets/,/pypars/,/pydat1/,/pydat2/,/pydat3/,/pysubs/,
 445 c--     +     /hyjets/,/parimp/,/hyjpar/,/hyipar/,/hyfpar/,/hypyin/
 446
 447 * generate 'njet' PYTHIA events and fill arrays for partons and hadrons
 448
 449 c       write(*,*)'in hyhard 0'
 450
 451       nshad=0
 452       noquen=0
 453       if(nhsel.eq.1.or.nhsel.eq.3) noquen=1
 454       if(njet.ge.1) then
 455        mdcy(pycomp(111),1)=0                     ! no pi0 decay
 456        mdcy(pycomp(310),1)=0                     ! no K_S0 decay
 457
 458        ifbp=0                                    ! fix impact parameter
 459        bfixp=RA*bgen
 460        Ap=AW                                     ! atomic weight
 461
 462        do ihard=1,njet
 463         mstp(111)=0
 464
 465 * generate type of NN sub-collision (pp, pn or nn)
 466         rand1=pyr(0)
 467         if(rand1.lt.rzta) then
 468          beam='p'
 469         else
 470          beam='n'
 471         end if
 472         rand2=pyr(0)
 473         if(rand2.lt.rzta) then
 474          targ='p'
 475         else
 476          targ='n'
 477         end if
 478         call pyinit('cms',beam,targ,ene)
 479 c        mstj(41)=0                           ! vacuum showering off
 480
 481
 482 c         write(*,*)'in hyhard 1'
 483
 484         call pyevnt                           ! generate hard scattering
 485
 486 c         write(*,*)'in hyhard 2'
 487
 488
 489 * PYQUEN: quenched jets if noquen=0 or non-quenched jets if noquen=1 and ishad=1
 490         if(ishad.eq.1.or.nhsel.eq.2.or.nhsel.eq.4)
 491      >   call pyquen(Ap,ifbp,bfixp)
 492
 493 c-ml    coordinate info if we need of it
 494
 495 c        write(*,*)'in hyhard pari(21)',pari(21)
 496
 497         Q=pari(21)
 498         x=r0*cos(psib1)         !fm
 499         y=r0*sin(psib1)  !fm
 500 c        write(*,*)'in hyhard ',x,y
 501
 502         if(Q.gt.0)then
 503         tau=1./Q ! 1/GeV
 504         tau=tau*0.197 !fm/c
 505         rm=0.0d0
 506         sig=1.0d0
 507
 508 c       write(*,*)'in hyhard tau',tau
 509
 510         etaLj=gaus(rm,sig) !fm
 511
 512 c       write(*,*)'in hyhard after gauss',tau
 513         z=tau*sinh(etaLj)
 514         t=tau*cosh(etaLj)
 515 c       write(*,*)'in hyhard z t',z,t
 516         endif
 517
 518 c       write(*,*)'x y z t', x,y,z,t,etaLj,Q
 519
 520 * treatment of "nuclear shadowing" (for Pb, Au, Pd or Ca beams only)
 521
 522
 523         if(ishad.eq.1) then
 524          kfh1=abs(k(3,2))
 525          kfh2=abs(k(4,2))
 526          xh1=pari(33)
 527          xh2=pari(34)
 528          Q2=pari(22)
 529          shad=shad1(kfh1,xh1,Q2,rb1)*shad1(kfh2,xh2,Q2,rb2)
 530          if(pyr(0).gt.shad) then
 531           nshad=nshad+1
 532           goto 53
 533          end if
 534         end if
 535
 536
 537 c       write(*,*)'in hyhard 4'
 538
 539         call pyexec                         ! hadronization done
 540
 541 c        write(*,*)'in hyhard 5'
 542
 543 c-ml
 544         call pyedit(2)                      ! remove partons & leave hadrons
 545
 546 c        write(*,*)'in hyhard 6'
 547
 548
 549 * fill array of final particles
 550         nu=nhj+n
 551         if(nu.gt.150000-np) then
 552          write(6,*) 'Warning, multiplicity too large! Cut hard part.'
 553          goto 52
 554         end if
 555         nhj=nu
 556
 557
 558         do i=nhj-n+1,nhj
 559          ip=i+n-nhj
 560 c           write(*,*)'in hyhard 7'
 561
 562          do j=1,5
 563           phj(i,j)=p(ip,j)
 564          end do
 565
 566           vhj(i,1)=x
 567           vhj(i,2)=y
 568           vhj(i,3)=z
 569           vhj(i,4)=t
 570
 571          do j=1,5
 572           khj(i,j)=k(ip,j)
 573          end do
 574          do j=3,5
 575           kk=khj(i,j)
 576           if(kk.gt.0) then
 577            khj(i,j)=kk+nhj-n
 578           end if
 579          end do
 580         end do
 581
 582
 583  53     continue
 584
 585 c       write(*,*)'in hyhard 8'
 586
 587        end do
 588  52    njet=ihard-1
 589
 590 c       write(*,*)'in hyhard 9', njet,nshad
 591
 592       end if
 593       njet=njet-nshad
 594
 595
 596 c       write(*,*)'in hyhard OUT'
 597
 598
 599       return
 600       end
 601 ****************************** END HYHARD **************************
 602 ********************************* HIPSEAR ***************************
 603       SUBROUTINE HIPSEAR (fmax,xmin)
 604 * find maximum and 'sufficient minimum' of differential inelasic AA cross
 605 * section as a function of impact paramater (xm, fm are outputs)
 606       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
 607       double precision npar0,nbco0
 608       external hsin
 609       common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
 610 c--      save /hyipar/
 611       xmin=(bmaxh-bminh)*RA
 612
 613 c      write(*,*)'bmaxh bminh',bmaxh,bminh
 614
 615       fmax=0.d0
 616       do j=1,1000
 617       x=bminh*RA+xmin*(j-1)/999.d0
 618       f=hsin(x)
 619        if(f.gt.fmax) then
 620         fmax=f
 621        endif
 622       end do
 623       return
 624       end
 625 ****************************** END HIPSEAR **************************
 626
 627 * differential inelastic AA cross section
 628       double precision function hsin(x)
 629       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
 630       external hftaa
 631       common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet
 632 c--      save /hyjpar/
 633       br=x
 634       hsin=br*(1.d0-dexp(-0.1d0*hftaa(br)*sigin))
 635       return
 636       end
 637
 638 * number of nucleons-participants at impact parameter b
 639       double precision function numpar(c)
 640       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
 641       external HFUNC1
 642       common /hynup1/ bp,x
 643       EPS=0.005d0
 644       A=0.d0
 645       B=6.28318d0
 646       H=0.01d0*(B-A)
 647       bp=c
 648       CALL SIMPA(A,B,H,EPS,1.d-8,HFUNC1,X,RES,AIH,AIABS)
 649       numpar=2.d0*RES
 650       return
 651       end
 652 *
 653       double precision function HFUNC1(x)
 654       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
 655       double precision npar0,nbco0
 656       external HFUNC2
 657       common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
 658       common /hynup1/ bp,xx
 659 c--      save /hyipar/
 660       if(init.eq.1) then
 661        Rl=Rpb
 662       else
 663        Rl=RA
 664       end if
 665       xx=x
 666       EPS=0.005d0
 667       A=0.d0
 668       B=3.d0*Rl
 669       H=0.01d0*(B-A)
 670       CALL SIMPB(A,B,H,EPS,1.d-8,HFUNC2,Y,RES,AIH,AIABS)
 671       HFUNC1=RES
 672       return
 673       end
 674 *
 675       double precision function HFUNC2(y)
 676       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
 677       double precision npar0,nbco0
 678       external hythik
 679       common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
 680       common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet
 681       common /hynup1/ bp,x
 682 c--      save /hyipar/,/hyjpar/
 683       r1=dsqrt(abs(y*y+bp*bp/4.d0+y*bp*dcos(x)))
 684       r2=dsqrt(abs(y*y+bp*bp/4.d0-y*bp*dcos(x)))
 685       s=1.d0-dexp(-0.1d0*sigin*hythik(r2))
 686       HFUNC2=y*hythik(r1)*s
 687       return
 688       end
 689
 690 * nuclear overlap function at impact parameter b
 691       double precision function hftaa(c)
 692       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
 693       common /hythic/ BAB(110),TAB(110),TAAB(110)
 694 c--      save /hythic/
 695       call parinv(c,BAB,TAAB,110,RES)
 696       hftaa=RES
 697       return
 698       end
 699 *
 700       double precision function HFUNC3(x)
 701       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
 702       double precision npar0,nbco0
 703       external HFUNC4
 704       common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
 705       common /hynup1/ bp,xx
 706 c--      save /hyipar/
 707       if(init.eq.1) then
 708        Rl=Rpb
 709       else
 710        Rl=RA
 711       end if
 712       xx=x
 713       EPS=0.005d0
 714       A=0.d0
 715       B=3.d0*Rl
 716       H=0.01d0*(B-A)
 717       CALL SIMPB(A,B,H,EPS,1.d-8,HFUNC4,Y,RES,AIH,AIABS)
 718       HFUNC3=RES
 719       return
 720       end
 721 *
 722       double precision function HFUNC4(y)
 723       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
 724       double precision npar0,nbco0
 725       external hythik
 726       common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
 727       common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet
 728       common /hynup1/ bp,x
 729 c--      save /hyipar/,/hyjpar/
 730       r1=dsqrt(abs(y*y+bp*bp/4.d0+y*bp*dcos(x)))
 731       r2=dsqrt(abs(y*y+bp*bp/4.d0-y*bp*dcos(x)))
 732       HFUNC4=y*hythik(r1)*hythik(r2)
 733       return
 734       end
 735
 736 * nuclear thickness function
 737        double precision function hythik(r)
 738        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
 739        common /hythic/ BAB(110),TAB(110),TAAB(110)
 740 c--       save /hythic/
 741        call parinv(r,BAB,TAB,110,RES)
 742        hythik=RES
 743        return
 744        end
 745
 746 * Wood-Saxon nucleon distrubution
 747        double precision function rhoaa(z)
 748        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
 749        double precision npar0,nbco0
 750        common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
 751        common /hygeom/ BC
 752 c--       save /hyipar/,/hygeom/
 753        if(init.eq.1) then
 754         Rl=Rpb
 755        else
 756         Rl=RA
 757        end if
 758        pi=3.14159d0
 759        df=0.54d0
 760        r=sqrt(bc*bc+z*z)
 761        rho0=3.d0/(4.d0*pi*Rl**3)/(1.d0+(pi*df/Rl)**2)
 762        rhoaa=rho0/(1.d0+dexp((r-Rl)/df))
 763        return
 764        end
 765
 766 * function to calculate nuclear shadowing factor (for Pb, Au, Pd or Ca beams)
 767       double precision function shad1(kfh,xbjh,Q2h,r)
 768       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
 769       double precision npar0,nbco0,nbcol,npart,npart0
 770       external ggshad
 771       common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
 772       common /hyfpar/ bgen,nbcol,npart,npart0,npyt,nhyd
 773       common /hyshad/ bbmin,bbmax,inuc
 774 c--      save /hyipar/,/hyfpar/,/hyshad/
 775       dimension res(2)
 776       kf=kfh
 777       xbj=xbjh
 778       Q2=Q2h
 779       bb=r
 780       inuc=0
 781       if(AW.gt.205.d0.and.AW.lt.209.d0) inuc=4         ! Pb-206, 207 or 208
 782       if(AW.gt.196.d0.and.AW.lt.198.d0) inuc=3         ! Au-197
 783       if(AW.gt.109.d0.and.AW.lt.111.d0) inuc=2         ! Pd-110
 784       if(AW.gt.39.d0.and.AW.lt.41.d0) inuc=1           ! Ca-40
 785       if(inuc.eq.0.and.ipr.eq.1) then
 786        write(6,*)
 787      > 'Warning! Shadowing is not foreseen for atomic weigth  A  ='
 788      > ,AW
 789        write(6,*)'******************************************************
 790      >************************'
 791        ipr=0
 792       end if
 793       if(inuc.ne.0) then
 794        xbj=max(5.d-5,xbj)
 795        xbj=min(0.95d0,xbj)
 796        Q2=max(4.d0,Q2)
 797        Q2=min(520.d0,Q2)
 798        bb=max(0.d0,bb)
 799        call ggshad(inuc,xbj,Q2,bb,res,ta)
 800        if(kf.eq.21) then
 801         shad1=res(1)
 802        elseif(kf.eq.1.or.kf.eq.2.or.kf.eq.3) then
 803         shad1=res(2)
 804        else
 805         shad1=1.d0
 806        end if
 807        else
 808         shad1=1.d0
 809       end if
 810       return
 811       end
 812
 813 ******************************************************************************
 814 * The part of the code which follows below, includes nuclear shadowing model *
 815 * and has been written by Konrad Tywoniuk (Oslo University, Norway)          *
 816 ******************************************************************************
 817 c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$
 818 c$$$
 819 c$$$
 820 c$$$  Shadowing from Glauber-Gribov theory for
 821 c$$$  gluons and light quarks (u,d,s and their antiquarks).
 822 c$$$  All Pomeron tree diagrams have been summed (Schwimmer model).
 823 c$$$  More details about the model in
 824 c$$$  K. Tywoniuk, I.C. Arsene, L. Bravina, A. Kaidalov and E.
 825 c$$$  Zabrodin, Phys. Lett. B 657 (2007) 170
 826 c$$$
 827 c$$$  We use FIT B from the H1 parameterization published in
 828 c$$$  A. Aktas et al. (H1 Collaboration), Eur. Phys. J C 48 (2006) 715
 829 c$$$  A. Aktas et al. (H1 Collaboration), Eur. Phys. J C 48 (2006) 749
 830 c$$$
 831 c$$$  Main routine is GGSHAD which provides the user with the
 832 c$$$  ratio of nuclear and nucleon parton distribution function
 833 c$$$  normalized by the atomic number for a given
 834 c$$$  INUCL:     nucleus (1=Ca,2=Pd,3=Au,4=Pb)
 835 c$$$  X:         Bjorken x
 836 c$$$  Q2:        q**2, scale squared
 837 c$$$  B:         impact parameter
 838 c$$$
 839 c$$$  Then      RES(1):    gluon shadowing
 840 c$$$            RES(2):    sea quark shadowing
 841 c$$$            TA:        nuclear profile function
 842 c$$$
 843 c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$
 844       SUBROUTINE GGSHAD(INUCL,X,Q2,B,RES,TAF)
 845       IMPLICIT NONE
 846       DOUBLE PRECISION TA(31,4),IMPAR(31),ANUCL(4)
 847       DOUBLE PRECISION XB(36),Q2V(13)
 848       DOUBLE PRECISION XMAX,XMAXX,Q2MIN,Q2MAX,BMAX
 849       DOUBLE PRECISION G(36,13,4),LQ(36,13,4)
 850       DOUBLE PRECISION TATMP(31),SHAD(2)
 851       DOUBLE PRECISION C(100),D(100),E(100)
 852       DOUBLE PRECISION X,Q2,B
 853       DOUBLE PRECISION RES(2)
 854       DOUBLE PRECISION TAF
 855       DOUBLE PRECISION SEVAL
 856       INTEGER INUCL,IREAD,TMAX
 857       INTEGER I,J,K,IK
 858
 859       PARAMETER(TMAX=31)
 860       PARAMETER(XMAX=0.1)
 861       PARAMETER(XMAXX=0.95)
 862       PARAMETER(Q2MIN=4.)
 863       PARAMETER(Q2MAX=520.)
 864       PARAMETER(BMAX=30.)
 865
 866       COMMON/GG07/XB,Q2V,G,LQ
 867
 868       DATA ANUCL/40.,110.,197.,206./
 869       DATA IMPAR
 870      >     /0.,0.5,1.,1.5,2.,2.5,3.,3.5,4.,4.5,5.,5.5,6.,6.5,7.,7.5,8.,
 871      >     8.5,9.,9.5,10.,10.5,11.,11.5,12.,12.5,13.,13.5,14.,14.5,15./
 872       DATA TA
 873      >     /0.0291662,0.0288636,0.0279354,0.0263044,0.0238341,0.0203565,
 874      >     0.0158333,0.0107343,0.00617758,0.00307859,0.00139738,
 875      >     0.000603839,0.000254843,0.000106329,4.40965e-05,1.82206e-05,
 876      >     7.50949e-06,3.08881e-06,1.26837e-06,5.20082e-07,2.12978e-07,
 877      >     8.71154e-08,3.55956e-08,1.45304e-08,5.92622e-09,2.41504e-09,
 878      >     9.83429e-10,4.00184e-10,1.62741e-10,6.61407e-11,2.68657e-11,
 879      >     0.0153487,0.0152775,0.0150614,0.0146926,0.0141561,0.0134268,
 880      >     0.0124651,0.0112115,0.00959258,0.00756907,0.00527514,
 881      >     0.00313482,0.00159852,0.000732479,0.000316628,0.000133116,
 882      >     5.52508e-05,2.27904e-05,9.3692e-06,3.84349e-06,1.57423e-06,
 883      >     6.43962e-07,2.63132e-07,1.07414e-07,4.38089e-08,1.78529e-08,
 884      >     7.2699e-09,2.95832e-09,1.20304e-09,4.8894e-10,1.98603e-10,
 885      >     0.0105299,0.010498,0.0104017,0.010239,0.010006,0.00969677,
 886      >     0.00930205,0.00880781,0.00819264,0.00742455,0.00646058,
 887      >     0.00526252,0.00385835,0.00243988,0.00131252,0.000621079,
 888      >     0.000272322,0.000114993,4.77312e-05,1.96569e-05,8.06377e-06,
 889      >     3.30066e-06,1.34904e-06,5.50751e-07,2.24633e-07,9.15434e-08,
 890      >     3.72777e-08,1.51694e-08,6.16885e-09,2.50714e-09,1.01838e-09,
 891      >     0.0102179,0.0101879,0.0100975,0.00994467,0.00972606,
 892      >     0.00943628,0.0090671,0.00860605,0.00803416,0.00732289,
 893      >     0.00643266,0.00532296,0.00400025,0.00261264,0.00144993,
 894      >     0.00070098,0.000310867,0.000131949,5.48887e-05,2.26247e-05,
 895      >     9.28457e-06,3.80091e-06,1.55358e-06,6.34273e-07,2.58701e-07,
 896      >     1.05427e-07,4.29316e-08,1.74701e-08,7.10447e-09,2.88739e-09,
 897      >     1.17283e-09/
 898       DATA Q2V
 899      >     /4.000,   6.000,   9.000,  13.500,  20.250,  30.375,  45.562,
 900      >     68.344, 102.516, 153.773, 230.660, 345.990, 518.985/
 901       DATA XB
 902      >     /0.99999998E-05,  0.13000000E-04,  0.16899999E-04,
 903      >     0.21970000E-04,  0.28561000E-04,  0.37129299E-04,
 904      >     0.48268099E-04,  0.62748499E-04,  0.81573096E-04,
 905      >     0.10604500E-03,  0.13785800E-03,  0.17921600E-03,
 906      >     0.23298099E-03,  0.30287501E-03,  0.39373801E-03,
 907      >     0.51185902E-03,  0.66541700E-03,  0.86504198E-03,
 908      >     0.11245500E-02,  0.14619200E-02,  0.19005000E-02,
 909      >     0.24706500E-02,  0.32118401E-02,  0.41753901E-02,
 910      >     0.54280101E-02,  0.70564100E-02,  0.91733299E-02,
 911      >     0.11925300E-01,  0.15502900E-01,  0.20153800E-01,
 912      >     0.26200000E-01,  0.34059901E-01,  0.44277899E-01,
 913      >     0.57561301E-01,  0.74829698E-01,  0.97278602E-01/
 914       DATA G
 915      >     /2.542140,2.457280,2.372200,2.283920,2.195320,2.105540,
 916      >     2.014040,1.921480,1.827600,1.733930,1.639990,1.546170,
 917      >     1.453220,1.361040,1.268750,1.178500,1.088730,1.000740,
 918      >     0.914606,0.830842,0.749315,0.671419,0.596883,0.526173,
 919      >     0.460067,0.398489,0.341895,0.289998,0.242628,0.199405,
 920      >     0.159625,0.122521,0.087252,0.053699,0.022804,0.000414,
 921      >     1.825600,1.768440,1.712550,1.654660,1.595830,1.533980,
 922      >     1.472860,1.409910,1.346060,1.282030,1.217520,1.153050,
 923      >     1.089030,1.025930,0.962709,0.900094,0.837787,0.776733,
 924      >     0.716072,0.656570,0.598295,0.541607,0.486661,0.433846,
 925      >     0.383421,0.335585,0.290685,0.248542,0.209305,0.172614,
 926      >     0.138166,0.105552,0.074246,0.044592,0.017935,0.000291,
 927      >     1.445600,1.405920,1.364150,1.322050,1.276900,1.232020,
 928      >     1.184670,1.136620,1.086950,1.037370,0.987225,0.937488,
 929      >     0.888244,0.839061,0.790370,0.741798,0.693426,0.645511,
 930      >     0.598298,0.551516,0.505500,0.460400,0.416326,0.373488,
 931      >     0.332281,0.292614,0.254859,0.218944,0.184976,0.152759,
 932      >     0.122077,0.092747,0.064477,0.037873,0.014527,0.000216,
 933      >     1.209700,1.179300,1.146560,1.112700,1.077790,1.041400,
 934      >     1.003150,0.964598,0.923825,0.882423,0.841304,0.800381,
 935      >     0.759486,0.718947,0.678587,0.638579,0.598671,0.558978,
 936      >     0.519524,0.480617,0.442116,0.404231,0.367057,0.330599,
 937      >     0.295291,0.261100,0.228220,0.196605,0.166367,0.137331,
 938      >     0.109496,0.082677,0.056835,0.032696,0.012013,0.000165,
 939      >     1.049390,1.024810,0.998839,0.971705,0.942633,0.912223,
 940      >     0.880504,0.847236,0.812716,0.777477,0.741670,0.706578,
 941      >     0.671021,0.636183,0.601384,0.566923,0.532196,0.497956,
 942      >     0.463884,0.430069,0.396612,0.363536,0.330915,0.298900,
 943      >     0.267627,0.237155,0.207697,0.179169,0.151566,0.124961,
 944      >     0.099266,0.074436,0.050544,0.028462,0.010013,0.000127,
 945      >     0.943086,0.922388,0.900700,0.877793,0.853134,0.826808,
 946      >     0.798937,0.769284,0.738709,0.706829,0.674818,0.643041,
 947      >     0.611250,0.580168,0.548982,0.517908,0.486720,0.456021,
 948      >     0.425312,0.394899,0.364616,0.334682,0.305173,0.276042,
 949      >     0.247465,0.219584,0.192474,0.166033,0.140391,0.115481,
 950      >     0.091347,0.068001,0.045636,0.025180,0.008510,0.000101,
 951      >     0.858923,0.841883,0.823827,0.803952,0.782124,0.759447,
 952      >     0.734659,0.708141,0.680281,0.651690,0.622372,0.593451,
 953      >     0.564610,0.536045,0.507521,0.479248,0.450845,0.422731,
 954      >     0.394687,0.366774,0.339095,0.311580,0.284350,0.257527,
 955      >     0.231101,0.205195,0.179929,0.155178,0.131077,0.107555,
 956      >     0.084712,0.062629,0.041557,0.022497,0.007326,0.000082,
 957      >     0.790086,0.776082,0.760408,0.743282,0.724441,0.703812,
 958      >     0.681673,0.657893,0.632354,0.605918,0.578976,0.552364,
 959      >     0.525820,0.499665,0.473306,0.447125,0.421027,0.394926,
 960      >     0.369017,0.343250,0.317596,0.292061,0.266805,0.241770,
 961      >     0.217121,0.192830,0.169123,0.145809,0.122984,0.100685,
 962      >     0.078967,0.057983,0.038065,0.020236,0.006364,0.000067,
 963      >     0.732957,0.720881,0.707633,0.692613,0.676200,0.657914,
 964      >     0.637626,0.615846,0.592373,0.567813,0.542877,0.518036,
 965      >     0.493386,0.468899,0.444596,0.420291,0.395908,0.371641,
 966      >     0.347423,0.323358,0.299396,0.275511,0.251830,0.228313,
 967      >     0.205146,0.182274,0.159813,0.137693,0.115983,0.094703,
 968      >     0.073974,0.053979,0.035068,0.018327,0.005580,0.000056,
 969      >     0.684126,0.674389,0.662432,0.649473,0.634783,0.618458,
 970      >     0.600173,0.580108,0.558464,0.535421,0.512051,0.488760,
 971      >     0.465717,0.442837,0.419997,0.397309,0.374393,0.351598,
 972      >     0.328943,0.306305,0.283685,0.261210,0.238875,0.216696,
 973      >     0.194737,0.173035,0.151692,0.130607,0.109864,0.089483,
 974      >     0.069628,0.050489,0.032492,0.016714,0.004940,0.000047,
 975      >     0.642684,0.634012,0.624311,0.612687,0.599751,0.584885,
 976      >     0.568165,0.549716,0.529352,0.507582,0.485673,0.463844,
 977      >     0.442031,0.420480,0.398892,0.377446,0.355834,0.334383,
 978      >     0.312903,0.291520,0.270128,0.248790,0.227632,0.206518,
 979      >     0.185624,0.164939,0.144527,0.124348,0.104415,0.084843,
 980      >     0.065754,0.047403,0.030217,0.015305,0.004396,0.000040,
 981      >     0.606784,0.599501,0.591216,0.581053,0.569262,0.555644,
 982      >     0.540298,0.522979,0.503888,0.483449,0.462725,0.441941,
 983      >     0.421308,0.400965,0.380515,0.360166,0.339667,0.319306,
 984      >     0.298896,0.278531,0.258224,0.237932,0.217694,0.197574,
 985      >     0.177567,0.157731,0.138168,0.118766,0.099587,0.080708,
 986      >     0.062310,0.044656,0.028210,0.014081,0.003935,0.000035,
 987      >     0.575250,0.569055,0.561834,0.552756,0.542358,0.530173,
 988      >     0.515880,0.499697,0.481592,0.462197,0.442485,0.422845,
 989      >     0.403137,0.383704,0.364332,0.344943,0.325467,0.306020,
 990      >     0.286496,0.267106,0.247644,0.228251,0.208881,0.189557,
 991      >     0.170377,0.151316,0.132475,0.113748,0.095226,0.076992,
 992      >     0.059214,0.042195,0.026424,0.013004,0.003541,0.000030,
 993      >     2.510860,2.424890,2.338800,2.251430,2.161030,2.070080,
 994      >     1.976500,1.884380,1.789520,1.695430,1.600900,1.507110,
 995      >     1.413980,1.321430,1.229430,1.138950,1.049380,0.961914,
 996      >     0.875863,0.792403,0.711728,0.634220,0.560088,0.490193,
 997      >     0.424574,0.363647,0.307681,0.256400,0.209654,0.167105,
 998      >     0.128273,0.092911,0.060892,0.033178,0.011804,0.000172,
 999      >     1.797890,1.740890,1.683730,1.625330,1.564220,1.502940,
1000      >     1.439440,1.375560,1.311540,1.247060,1.182260,1.117820,
1001      >     1.053940,0.990598,0.927193,0.864982,0.802868,0.741596,
1002      >     0.681232,0.621931,0.563990,0.507817,0.452982,0.400714,
1003      >     0.350809,0.303514,0.259139,0.217505,0.178849,0.142987,
1004      >     0.109689,0.079022,0.051121,0.027212,0.009208,0.000121,
1005      >     1.421860,1.380980,1.338550,1.294150,1.249120,1.201850,
1006      >     1.153940,1.105280,1.054910,1.004680,0.954677,0.905086,
1007      >     0.855141,0.806467,0.757283,0.708817,0.660443,0.613233,
1008      >     0.565935,0.519307,0.473770,0.428741,0.385002,0.342676,
1009      >     0.301797,0.262591,0.225442,0.189974,0.156683,0.125290,
1010      >     0.095921,0.068630,0.043899,0.022874,0.007411,0.000090,
1011      >     1.187690,1.155530,1.122280,1.087030,1.051020,1.013270,
1012      >     0.973852,0.933953,0.892905,0.851713,0.810455,0.769192,
1013      >     0.728379,0.687798,0.647423,0.607517,0.567475,0.528048,
1014      >     0.488850,0.450119,0.411897,0.374353,0.337466,0.301396,
1015      >     0.266522,0.232713,0.200387,0.169339,0.139725,0.111661,
1016      >     0.085284,0.060612,0.038324,0.019572,0.006096,0.000069,
1017      >     1.028240,1.002610,0.975714,0.947100,0.916605,0.885287,
1018      >     0.852214,0.817988,0.782800,0.747390,0.711939,0.676325,
1019      >     0.641102,0.606200,0.571466,0.536994,0.502546,0.468438,
1020      >     0.434547,0.400992,0.367663,0.334910,0.302578,0.270890,
1021      >     0.240102,0.210041,0.181153,0.153175,0.126406,0.100826,
1022      >     0.076605,0.054074,0.033774,0.016896,0.005055,0.000053,
1023      >     0.923239,0.901357,0.878307,0.853928,0.827994,0.800448,
1024      >     0.770833,0.740864,0.709208,0.677436,0.645505,0.613815,
1025      >     0.582211,0.551121,0.519863,0.488847,0.457919,0.427291,
1026      >     0.396728,0.366593,0.336576,0.306991,0.277689,0.248974,
1027      >     0.220914,0.193485,0.166892,0.141077,0.116297,0.092513,
1028      >     0.069960,0.049024,0.030239,0.014840,0.004277,0.000042,
1029      >     0.840109,0.821572,0.802221,0.780737,0.758007,0.733570,
1030      >     0.707391,0.680060,0.651741,0.622598,0.593525,0.564790,
1031      >     0.535923,0.507607,0.479180,0.450770,0.422703,0.394808,
1032      >     0.366922,0.339201,0.311718,0.284555,0.257727,0.231240,
1033      >     0.205291,0.179821,0.155134,0.131145,0.107877,0.085551,
1034      >     0.064423,0.044815,0.027334,0.013173,0.003668,0.000034,
1035      >     0.771898,0.756403,0.739454,0.720917,0.700462,0.678667,
1036      >     0.655360,0.630295,0.604073,0.577412,0.550894,0.524237,
1037      >     0.497613,0.471552,0.445325,0.419311,0.393354,0.367582,
1038      >     0.341816,0.316292,0.290960,0.265746,0.240859,0.216190,
1039      >     0.191974,0.168255,0.145100,0.122534,0.100715,0.079595,
1040      >     0.059666,0.041209,0.024869,0.011778,0.003176,0.000028,
1041      >     0.715432,0.701697,0.687335,0.671043,0.652774,0.633196,
1042      >     0.611618,0.588703,0.564555,0.539941,0.515046,0.490273,
1043      >     0.465666,0.441417,0.417189,0.393083,0.368842,0.344918,
1044      >     0.320853,0.296990,0.273346,0.249853,0.226554,0.203403,
1045      >     0.180727,0.158303,0.136484,0.115160,0.094448,0.074509,
1046      >     0.055551,0.038122,0.022769,0.010610,0.002775,0.000023,
1047      >     0.666960,0.655627,0.643169,0.628497,0.612317,0.594375,
1048      >     0.574670,0.553372,0.531024,0.507805,0.484742,0.461630,
1049      >     0.438595,0.415762,0.393171,0.370479,0.347803,0.325380,
1050      >     0.302957,0.280541,0.258199,0.236074,0.214141,0.192375,
1051      >     0.170880,0.149710,0.129011,0.108739,0.089084,0.069999,
1052      >     0.052020,0.035461,0.020980,0.009629,0.002450,0.000020,
1053      >     0.626089,0.616427,0.605000,0.592191,0.577470,0.561221,
1054      >     0.543098,0.523316,0.502350,0.480560,0.458749,0.436948,
1055      >     0.415413,0.393933,0.372648,0.351193,0.329853,0.308661,
1056      >     0.287323,0.266214,0.245090,0.224191,0.203371,0.182736,
1057      >     0.162324,0.142202,0.122478,0.103078,0.084237,0.066064,
1058      >     0.048902,0.033114,0.019409,0.008778,0.002175,0.000017,
1059      >     0.590664,0.582220,0.572479,0.560768,0.547507,0.532490,
1060      >     0.515657,0.497052,0.477369,0.456893,0.436199,0.415686,
1061      >     0.395192,0.374819,0.354542,0.334339,0.314120,0.293919,
1062      >     0.273897,0.253708,0.233747,0.213772,0.193977,0.174220,
1063      >     0.154737,0.135524,0.116625,0.098102,0.080014,0.062548,
1064      >     0.046083,0.031041,0.018032,0.008042,0.001942,0.000014,
1065      >     0.559453,0.552473,0.543647,0.533124,0.521070,0.507229,
1066      >     0.491452,0.474143,0.455506,0.435945,0.416256,0.396977,
1067      >     0.377326,0.358052,0.338792,0.319598,0.300262,0.281122,
1068      >     0.261883,0.242706,0.223579,0.204566,0.185577,0.166710,
1069      >     0.148088,0.129570,0.111434,0.093593,0.076187,0.059415,
1070      >     0.043600,0.029214,0.016816,0.007398,0.001744,0.000012,
1071      >     2.494430,2.407750,2.320680,2.230570,2.139410,2.045820,
1072      >     1.954300,1.859160,1.764370,1.669990,1.574490,1.481120,
1073      >     1.386830,1.294720,1.202200,1.111670,1.021960,0.934479,
1074      >     0.848525,0.765125,0.684539,0.607178,0.533502,0.463759,
1075      >     0.398700,0.338188,0.282593,0.231728,0.185473,0.143600,
1076      >     0.105860,0.072262,0.043334,0.020543,0.005838,0.000062,
1077      >     1.781830,1.724140,1.666660,1.605560,1.544230,1.481280,
1078      >     1.418440,1.353130,1.288210,1.223120,1.158330,1.093950,
1079      >     1.029490,0.965787,0.902413,0.840091,0.777784,0.716620,
1080      >     0.656180,0.597156,0.539338,0.483075,0.428748,0.376710,
1081      >     0.327154,0.280125,0.236105,0.194961,0.156820,0.121629,
1082      >     0.089495,0.060685,0.035892,0.016616,0.004508,0.000044,
1083      >     1.407480,1.365110,1.320730,1.276090,1.228930,1.181910,
1084      >     1.132930,1.082570,1.032770,0.982143,0.932084,0.882066,
1085      >     0.832465,0.783250,0.734272,0.685821,0.637505,0.589736,
1086      >     0.542652,0.496248,0.450545,0.405906,0.362451,0.320187,
1087      >     0.279650,0.240757,0.203901,0.169026,0.136216,0.105650,
1088      >     0.077490,0.052176,0.030469,0.013808,0.003599,0.000032,
1089      >     1.173640,1.140320,1.105880,1.069540,1.032320,0.993897,
1090      >     0.953913,0.913021,0.871617,0.829958,0.788517,0.747238,
1091      >     0.706258,0.666094,0.625534,0.585422,0.545512,0.506014,
1092      >     0.466972,0.428374,0.390208,0.352773,0.316057,0.280282,
1093      >     0.245648,0.212195,0.180160,0.149604,0.120645,0.093445,
1094      >     0.068267,0.045631,0.026330,0.011698,0.002940,0.000025,
1095      >     1.014900,0.987942,0.959684,0.930270,0.899044,0.866678,
1096      >     0.832642,0.797620,0.762091,0.726232,0.690987,0.655077,
1097      >     0.620091,0.585191,0.550543,0.515986,0.481432,0.447342,
1098      >     0.413584,0.380034,0.346933,0.314295,0.282265,0.250834,
1099      >     0.220232,0.190535,0.162003,0.134530,0.108400,0.083726,
1100      >     0.060878,0.040604,0.022978,0.010002,0.002422,0.000019,
1101      >     0.909560,0.887017,0.862965,0.837572,0.810211,0.781875,
1102      >     0.751746,0.720666,0.688703,0.656776,0.624847,0.593122,
1103      >     0.561398,0.530391,0.499138,0.468280,0.437229,0.406831,
1104      >     0.376391,0.346350,0.316458,0.286984,0.257979,0.229524,
1105      >     0.201720,0.174606,0.148473,0.123248,0.099111,0.076308,
1106      >     0.055184,0.036277,0.020402,0.008711,0.002038,0.000015,
1107      >     0.827323,0.808240,0.787194,0.764785,0.741102,0.715710,
1108      >     0.688724,0.660472,0.631563,0.602307,0.573344,0.544275,
1109      >     0.515643,0.487203,0.458901,0.430660,0.402505,0.374606,
1110      >     0.346855,0.319408,0.292040,0.265141,0.238553,0.212326,
1111      >     0.186666,0.161641,0.137398,0.113940,0.091458,0.070180,
1112      >     0.050481,0.032923,0.018295,0.007672,0.001739,0.000012,
1113      >     0.759389,0.742705,0.724897,0.705170,0.684041,0.661178,
1114      >     0.636646,0.610897,0.584460,0.557556,0.530715,0.504308,
1115      >     0.477779,0.451618,0.425597,0.399640,0.373702,0.348041,
1116      >     0.322356,0.297020,0.271779,0.246872,0.222152,0.197859,
1117      >     0.173964,0.150642,0.127976,0.105994,0.084903,0.064934,
1118      >     0.046463,0.030071,0.016527,0.006814,0.001499,0.000010,
1119      >     0.703019,0.688814,0.673059,0.655722,0.636535,0.615667,
1120      >     0.593245,0.569605,0.545086,0.520194,0.495411,0.470863,
1121      >     0.446201,0.421921,0.397793,0.373652,0.349528,0.325620,
1122      >     0.301845,0.278188,0.254649,0.231349,0.208318,0.185538,
1123      >     0.163153,0.141239,0.119901,0.099172,0.079281,0.060422,
1124      >     0.043027,0.027647,0.015032,0.006099,0.001305,0.000008,
1125      >     0.655549,0.643053,0.628966,0.613794,0.596440,0.577313,
1126      >     0.556691,0.534653,0.511805,0.488542,0.465448,0.442353,
1127      >     0.419418,0.396807,0.374153,0.351526,0.328993,0.306561,
1128      >     0.284259,0.262079,0.240000,0.218086,0.196361,0.174927,
1129      >     0.153790,0.133102,0.112916,0.093252,0.074414,0.056530,
1130      >     0.040069,0.025572,0.013769,0.005504,0.001148,0.000007,
1131      >     0.614886,0.603957,0.591709,0.577655,0.561834,0.544400,
1132      >     0.525193,0.504891,0.483435,0.461525,0.439729,0.418166,
1133      >     0.396478,0.375064,0.353725,0.332565,0.311308,0.290188,
1134      >     0.269125,0.248171,0.227323,0.206600,0.186070,0.165736,
1135      >     0.145686,0.125993,0.106807,0.088112,0.070133,0.053112,
1136      >     0.037464,0.023751,0.012665,0.004991,0.001015,0.000006,
1137      >     0.579569,0.570010,0.559085,0.546465,0.532125,0.516026,
1138      >     0.498089,0.478853,0.458681,0.438173,0.417404,0.396958,
1139      >     0.376499,0.356301,0.336115,0.316014,0.295930,0.275886,
1140      >     0.255871,0.236010,0.216223,0.196524,0.176991,0.157626,
1141      >     0.138515,0.119742,0.101404,0.083545,0.066339,0.050078,
1142      >     0.035167,0.022152,0.011703,0.004550,0.000904,0.000005,
1143      >     0.548753,0.540580,0.530520,0.519276,0.506004,0.491061,
1144      >     0.474436,0.456132,0.436996,0.417488,0.397755,0.378421,
1145      >     0.359085,0.339851,0.320700,0.301518,0.282357,0.263296,
1146      >     0.244273,0.225311,0.206421,0.187624,0.168959,0.150424,
1147      >     0.132145,0.114172,0.096557,0.079446,0.062967,0.047378,
1148      >     0.033126,0.020736,0.010987,0.004166,0.000809,0.000004,
1149      >     2.486670,2.401390,2.313150,2.223100,2.132210,2.040620,
1150      >     1.946240,1.852960,1.758660,1.663790,1.569300,1.475400,
1151      >     1.381630,1.288980,1.197380,1.106900,1.017420,0.930068,
1152      >     0.844430,0.760955,0.680549,0.603582,0.530075,0.460570,
1153      >     0.395641,0.335338,0.279919,0.229212,0.183101,0.141385,
1154      >     0.103826,0.070754,0.041898,0.019580,0.005430,0.000056,
1155      >     1.776150,1.719900,1.661320,1.600440,1.539110,1.477010,
1156      >     1.412790,1.348550,1.283220,1.218900,1.153830,1.088770,
1157      >     1.025080,0.961812,0.898582,0.836002,0.773940,0.712963,
1158      >     0.652650,0.593610,0.536099,0.480001,0.425916,0.373940,
1159      >     0.324454,0.277574,0.233716,0.192654,0.154624,0.119622,
1160      >     0.087676,0.059112,0.034659,0.015816,0.004189,0.000039,
1161      >     1.403340,1.361000,1.316930,1.271610,1.225380,1.177160,
1162      >     1.128240,1.078580,1.028640,0.977947,0.928063,0.878521,
1163      >     0.828550,0.779706,0.730569,0.682199,0.633948,0.586522,
1164      >     0.539401,0.493171,0.447639,0.403145,0.359673,0.317729,
1165      >     0.277222,0.238454,0.201716,0.166913,0.134251,0.103840,
1166      >     0.075864,0.050772,0.029388,0.013129,0.003341,0.000029,
1167      >     1.169950,1.136380,1.102100,1.065820,1.028190,0.989956,
1168      >     0.950189,0.908974,0.867893,0.826368,0.784940,0.743815,
1169      >     0.702872,0.662550,0.622187,0.582284,0.542431,0.503183,
1170      >     0.464112,0.425583,0.387521,0.350185,0.313576,0.277943,
1171      >     0.243398,0.210036,0.178125,0.147636,0.118821,0.091774,
1172      >     0.066794,0.044372,0.025371,0.011112,0.002727,0.000022,
1173      >     1.011370,0.984724,0.956435,0.926814,0.895443,0.862964,
1174      >     0.829274,0.794479,0.758875,0.722999,0.687335,0.652160,
1175      >     0.616882,0.582106,0.547373,0.512881,0.478596,0.444627,
1176      >     0.410813,0.377446,0.344502,0.311951,0.279953,0.248594,
1177      >     0.218124,0.188529,0.160091,0.132732,0.106696,0.082174,
1178      >     0.059499,0.039204,0.022122,0.009492,0.002245,0.000017,
1179      >     0.906689,0.883943,0.860162,0.834079,0.807067,0.778589,
1180      >     0.748446,0.717460,0.685716,0.653595,0.621726,0.590095,
1181      >     0.558631,0.527321,0.496368,0.465511,0.434680,0.404249,
1182      >     0.373821,0.343908,0.314123,0.284707,0.255819,0.227414,
1183      >     0.199673,0.172714,0.146637,0.121510,0.097519,0.074849,
1184      >     0.053902,0.035215,0.019623,0.008259,0.001888,0.000014,
1185      >     0.824084,0.805081,0.784198,0.761980,0.737877,0.712401,
1186      >     0.685185,0.657419,0.628573,0.599438,0.570456,0.541554,
1187      >     0.512854,0.484478,0.456174,0.428019,0.399930,0.372215,
1188      >     0.344523,0.317092,0.289826,0.262998,0.236396,0.210320,
1189      >     0.184715,0.159801,0.135667,0.112306,0.089938,0.068791,
1190      >     0.049274,0.031936,0.017585,0.007270,0.001610,0.000011,
1191      >     0.756740,0.739909,0.722048,0.702184,0.681164,0.658293,
1192      >     0.633526,0.608140,0.581469,0.554572,0.527904,0.501517,
1193      >     0.475145,0.448970,0.422970,0.397121,0.371248,0.345564,
1194      >     0.320059,0.294788,0.269623,0.244755,0.220148,0.195923,
1195      >     0.172134,0.148871,0.126312,0.104417,0.083450,0.063624,
1196      >     0.045337,0.029156,0.015873,0.006450,0.001387,0.000009,
1197      >     0.700680,0.686067,0.670562,0.652740,0.633579,0.612848,
1198      >     0.590476,0.566706,0.542207,0.517556,0.492696,0.468105,
1199      >     0.443562,0.419350,0.395237,0.371204,0.347185,0.323339,
1200      >     0.299543,0.276011,0.252614,0.229355,0.206336,0.183677,
1201      >     0.161386,0.139527,0.118283,0.097686,0.077903,0.059183,
1202      >     0.041956,0.026789,0.014429,0.005771,0.001207,0.000008,
1203      >     0.652991,0.640812,0.626563,0.610980,0.593505,0.574586,
1204      >     0.553939,0.531965,0.509196,0.485938,0.462796,0.439793,
1205      >     0.416876,0.394269,0.371716,0.349126,0.326720,0.304355,
1206      >     0.282067,0.259966,0.237971,0.216140,0.194490,0.173138,
1207      >     0.152089,0.131483,0.111379,0.091839,0.073081,0.055347,
1208      >     0.039055,0.024764,0.013208,0.005204,0.001061,0.000006,
1209      >     0.612493,0.601527,0.589260,0.574987,0.559293,0.541793,
1210      >     0.522714,0.502071,0.480580,0.459021,0.437202,0.415483,
1211      >     0.394052,0.372703,0.351511,0.330271,0.309080,0.288012,
1212      >     0.267062,0.246156,0.225353,0.204703,0.184235,0.163954,
1213      >     0.144013,0.124421,0.105317,0.086709,0.068861,0.051973,
1214      >     0.036498,0.022991,0.012143,0.004716,0.000938,0.000005,
1215      >     0.577224,0.567718,0.556768,0.543992,0.529614,0.513264,
1216      >     0.495593,0.476332,0.456118,0.435643,0.414958,0.394475,
1217      >     0.374159,0.353948,0.333917,0.313860,0.293761,0.273842,
1218      >     0.253865,0.234066,0.214285,0.194711,0.175202,0.155932,
1219      >     0.136903,0.118198,0.099962,0.082192,0.065107,0.048988,
1220      >     0.034247,0.021430,0.011215,0.004297,0.000835,0.000005,
1221      >     0.546580,0.538342,0.528410,0.516932,0.503593,0.488594,
1222      >     0.471844,0.453700,0.434588,0.415081,0.395496,0.376054,
1223      >     0.356702,0.337597,0.318488,0.299383,0.280225,0.261271,
1224      >     0.242274,0.223404,0.204576,0.185839,0.167243,0.148787,
1225      >     0.130586,0.112676,0.095184,0.078140,0.061781,0.046329,
1226      >     0.032244,0.020050,0.010400,0.003932,0.000747,0.000004/
1227       DATA LQ
1228      >     /0.514104,0.515832,0.517091,0.517536,0.517832,0.517221,
1229      >     0.516574,0.515224,0.513287,0.510858,0.508000,0.504917,
1230      >     0.501246,0.497403,0.492983,0.488152,0.482358,0.476122,
1231      >     0.468649,0.460339,0.450645,0.439812,0.427166,0.412836,
1232      >     0.396439,0.377753,0.356388,0.331566,0.302870,0.269780,
1233      >     0.231439,0.187668,0.138736,0.086850,0.036211,0.000604,
1234      >     0.579741,0.577476,0.575013,0.571386,0.567464,0.562828,
1235      >     0.557414,0.551410,0.544332,0.536891,0.528953,0.521049,
1236      >     0.512612,0.503967,0.494945,0.485619,0.475531,0.464958,
1237      >     0.453985,0.442176,0.429817,0.416392,0.402026,0.386267,
1238      >     0.369251,0.350413,0.329631,0.306106,0.279316,0.248539,
1239      >     0.213030,0.172492,0.127054,0.078865,0.032177,0.000513,
1240      >     0.617205,0.612889,0.607448,0.601280,0.593987,0.585725,
1241      >     0.577239,0.567645,0.557436,0.546468,0.535124,0.523690,
1242      >     0.511913,0.500192,0.488271,0.476132,0.463482,0.450653,
1243      >     0.437460,0.423920,0.409767,0.395178,0.379804,0.363620,
1244      >     0.346492,0.327958,0.307882,0.285606,0.260283,0.231441,
1245      >     0.198228,0.160224,0.117620,0.072457,0.029003,0.000445,
1246      >     0.635918,0.629589,0.622092,0.614001,0.604751,0.594717,
1247      >     0.583683,0.572041,0.559188,0.546115,0.532523,0.518899,
1248      >     0.505201,0.491782,0.477675,0.463752,0.449738,0.435623,
1249      >     0.421062,0.406484,0.391587,0.376279,0.360631,0.344349,
1250      >     0.327299,0.309280,0.289970,0.268596,0.244665,0.217365,
1251      >     0.185991,0.150104,0.109826,0.067166,0.026423,0.000392,
1252      >     0.642211,0.634584,0.625818,0.616705,0.605951,0.594718,
1253      >     0.582351,0.569037,0.554986,0.540276,0.525303,0.510521,
1254      >     0.495596,0.480637,0.465827,0.451058,0.435872,0.421066,
1255      >     0.405880,0.390796,0.375397,0.359903,0.344185,0.327987,
1256      >     0.311213,0.293614,0.274949,0.254454,0.231575,0.205516,
1257      >     0.175652,0.141477,0.103138,0.062618,0.024219,0.000348,
1258      >     0.639180,0.631092,0.621893,0.612137,0.600935,0.589006,
1259      >     0.576122,0.561754,0.547048,0.531226,0.515512,0.499878,
1260      >     0.484089,0.468662,0.453258,0.437735,0.422292,0.407026,
1261      >     0.391544,0.376258,0.360735,0.345208,0.329470,0.313579,
1262      >     0.297158,0.280066,0.261996,0.242200,0.220305,0.195335,
1263      >     0.166765,0.134023,0.097367,0.058717,0.022356,0.000312,
1264      >     0.631156,0.622788,0.613580,0.603527,0.592277,0.579779,
1265      >     0.566484,0.551813,0.536423,0.520508,0.504251,0.488021,
1266      >     0.472093,0.456329,0.440551,0.424908,0.409169,0.393745,
1267      >     0.378229,0.362803,0.347473,0.331992,0.316442,0.300771,
1268      >     0.284762,0.268174,0.250682,0.231645,0.210500,0.186498,
1269      >     0.158976,0.127569,0.092372,0.055348,0.020770,0.000282,
1270      >     0.619698,0.611519,0.602336,0.591943,0.580814,0.568514,
1271      >     0.554982,0.540509,0.524970,0.508623,0.492412,0.476046,
1272      >     0.459935,0.444036,0.428166,0.412462,0.396815,0.381305,
1273      >     0.365892,0.350546,0.335395,0.320121,0.304877,0.289437,
1274      >     0.273841,0.257707,0.240774,0.222305,0.201894,0.178726,
1275      >     0.152209,0.121898,0.087984,0.052399,0.019398,0.000257,
1276      >     0.606343,0.598499,0.589365,0.579821,0.568434,0.556327,
1277      >     0.542872,0.528454,0.512900,0.496651,0.480390,0.464149,
1278      >     0.447972,0.432310,0.416261,0.400748,0.385116,0.369727,
1279      >     0.354542,0.339434,0.324280,0.309328,0.294433,0.279362,
1280      >     0.264086,0.248363,0.231853,0.214025,0.194261,0.171833,
1281      >     0.146197,0.116876,0.084113,0.049805,0.018206,0.000236,
1282      >     0.591896,0.584256,0.575797,0.566248,0.555494,0.543797,
1283      >     0.530664,0.516396,0.501105,0.484951,0.468541,0.452368,
1284      >     0.436358,0.420764,0.405075,0.389617,0.374238,0.359020,
1285      >     0.343972,0.329043,0.314268,0.299607,0.284884,0.270198,
1286      >     0.255337,0.239949,0.223985,0.206644,0.187466,0.165701,
1287      >     0.140843,0.112412,0.080664,0.047512,0.017165,0.000218,
1288      >     0.576697,0.569585,0.561998,0.552972,0.542790,0.531174,
1289      >     0.518335,0.504287,0.489160,0.473285,0.457139,0.441230,
1290      >     0.425520,0.409903,0.394534,0.379177,0.364141,0.349046,
1291      >     0.334214,0.319602,0.305086,0.290690,0.276241,0.261890,
1292      >     0.247381,0.232411,0.216817,0.199969,0.181340,0.160182,
1293      >     0.135997,0.108375,0.077551,0.045446,0.016236,0.000202,
1294      >     0.561779,0.555491,0.548244,0.539719,0.529917,0.518731,
1295      >     0.506334,0.492791,0.477798,0.462049,0.446208,0.430534,
1296      >     0.414938,0.399619,0.384500,0.369572,0.354611,0.339912,
1297      >     0.325296,0.310898,0.296649,0.282465,0.268439,0.254375,
1298      >     0.240158,0.225578,0.210376,0.193910,0.175777,0.155185,
1299      >     0.131618,0.104727,0.074745,0.043586,0.015406,0.000188,
1300      >     0.547406,0.541522,0.534772,0.526637,0.517320,0.506665,
1301      >     0.494929,0.481541,0.466958,0.451581,0.435977,0.420382,
1302      >     0.405163,0.390039,0.375242,0.360562,0.345699,0.331352,
1303      >     0.316934,0.302875,0.288897,0.275075,0.261211,0.247494,
1304      >     0.233546,0.219335,0.204477,0.188465,0.170694,0.150613,
1305      >     0.127623,0.101395,0.072188,0.041898,0.014659,0.000176,
1306      >     0.510949,0.512461,0.513125,0.513775,0.513627,0.512953,
1307      >     0.511676,0.510153,0.508026,0.505586,0.502488,0.499514,
1308      >     0.495918,0.491777,0.487221,0.482318,0.476322,0.469865,
1309      >     0.462254,0.453539,0.443579,0.431957,0.418405,0.402714,
1310      >     0.384513,0.363211,0.338441,0.309525,0.276023,0.237418,
1311      >     0.194095,0.146991,0.098692,0.053921,0.018675,0.000251,
1312      >     0.573840,0.571539,0.568660,0.564636,0.560125,0.555429,
1313      >     0.549250,0.542654,0.535827,0.528339,0.520207,0.512007,
1314      >     0.503505,0.494854,0.485721,0.476289,0.466253,0.455649,
1315      >     0.444534,0.432471,0.419685,0.405709,0.390526,0.373656,
1316      >     0.355038,0.334208,0.310578,0.283443,0.252376,0.217005,
1317      >     0.177144,0.133976,0.089700,0.048645,0.016531,0.000213,
1318      >     0.610274,0.604837,0.598760,0.591967,0.584049,0.575825,
1319      >     0.566634,0.556589,0.545957,0.534847,0.523533,0.511998,
1320      >     0.500328,0.488702,0.476646,0.464532,0.451851,0.438900,
1321      >     0.425698,0.412022,0.397498,0.382650,0.366763,0.349574,
1322      >     0.330976,0.310680,0.288132,0.262614,0.233680,0.200673,
1323      >     0.163728,0.123628,0.082475,0.044455,0.014854,0.000185,
1324      >     0.627192,0.620322,0.611968,0.603362,0.593426,0.582683,
1325      >     0.571079,0.558798,0.545824,0.532347,0.519069,0.505318,
1326      >     0.491597,0.478201,0.464023,0.450314,0.436183,0.422097,
1327      >     0.407576,0.392963,0.377856,0.362363,0.346096,0.329071,
1328      >     0.310942,0.291304,0.269699,0.245618,0.218210,0.187279,
1329      >     0.152695,0.115069,0.076563,0.041003,0.013496,0.000163,
1330      >     0.631625,0.623836,0.614536,0.604467,0.593362,0.581189,
1331      >     0.568126,0.554473,0.540005,0.525342,0.510320,0.495556,
1332      >     0.480543,0.465856,0.450901,0.436047,0.421171,0.406166,
1333      >     0.391154,0.375993,0.360616,0.344868,0.328695,0.311826,
1334      >     0.294125,0.275057,0.254485,0.231491,0.205428,0.176164,
1335      >     0.143362,0.107930,0.071534,0.038049,0.012338,0.000145,
1336      >     0.628674,0.619289,0.609887,0.598892,0.587603,0.574461,
1337      >     0.560634,0.545957,0.530648,0.515019,0.499229,0.483499,
1338      >     0.467857,0.452587,0.437187,0.421817,0.406489,0.391076,
1339      >     0.375987,0.360550,0.345003,0.329399,0.313389,0.296812,
1340      >     0.279488,0.261191,0.241276,0.219209,0.194576,0.166551,
1341      >     0.135508,0.101696,0.067217,0.035540,0.011363,0.000130,
1342      >     0.619655,0.610842,0.600741,0.589506,0.577444,0.564098,
1343      >     0.550163,0.534861,0.519180,0.503427,0.486955,0.471011,
1344      >     0.455040,0.439305,0.423732,0.408132,0.392629,0.377169,
1345      >     0.361831,0.346498,0.331073,0.315724,0.299859,0.283638,
1346      >     0.266810,0.249063,0.229985,0.208716,0.185019,0.158309,
1347      >     0.128580,0.096355,0.063475,0.033376,0.010535,0.000117,
1348      >     0.607474,0.598863,0.589070,0.578071,0.565959,0.552353,
1349      >     0.538238,0.523141,0.507028,0.491211,0.474478,0.458296,
1350      >     0.442260,0.426619,0.410731,0.395250,0.379629,0.364356,
1351      >     0.348971,0.333841,0.318608,0.303265,0.287813,0.272025,
1352      >     0.255702,0.238406,0.219910,0.199549,0.176755,0.151092,
1353      >     0.122529,0.091666,0.060228,0.031496,0.009822,0.000107,
1354      >     0.593760,0.585285,0.576038,0.564716,0.553007,0.539720,
1355      >     0.525715,0.510544,0.494716,0.478059,0.461888,0.445882,
1356      >     0.429808,0.413922,0.398545,0.382937,0.367557,0.352300,
1357      >     0.337238,0.322197,0.307168,0.292182,0.277095,0.261674,
1358      >     0.245761,0.228939,0.211113,0.191442,0.169413,0.144749,
1359      >     0.117225,0.087581,0.057361,0.029844,0.009202,0.000098,
1360      >     0.579324,0.570883,0.561883,0.551297,0.539642,0.526840,
1361      >     0.512914,0.497903,0.481998,0.465838,0.449788,0.433693,
1362      >     0.417703,0.402286,0.386792,0.371437,0.356215,0.341283,
1363      >     0.326313,0.311637,0.296947,0.282140,0.267321,0.252270,
1364      >     0.236761,0.220580,0.203236,0.184230,0.162954,0.139137,
1365      >     0.112557,0.083990,0.054824,0.028391,0.008663,0.000090,
1366      >     0.564315,0.556562,0.547782,0.537905,0.526572,0.514011,
1367      >     0.500033,0.485378,0.470013,0.453915,0.438022,0.422153,
1368      >     0.406530,0.391112,0.375843,0.360792,0.345939,0.331089,
1369      >     0.316383,0.301913,0.287442,0.273075,0.258564,0.243883,
1370      >     0.228827,0.213013,0.196181,0.177712,0.157167,0.134036,
1371      >     0.108376,0.080670,0.052555,0.027085,0.008183,0.000084,
1372      >     0.549474,0.542540,0.534001,0.524222,0.513538,0.501429,
1373      >     0.487998,0.473665,0.458530,0.442760,0.426794,0.411187,
1374      >     0.395839,0.380807,0.365652,0.350917,0.336154,0.321565,
1375      >     0.307328,0.293072,0.278886,0.264801,0.250668,0.236339,
1376      >     0.221517,0.206147,0.189835,0.171902,0.151860,0.129450,
1377      >     0.104545,0.077707,0.050505,0.025913,0.007753,0.000078,
1378      >     0.535066,0.528539,0.520546,0.511350,0.500893,0.489256,
1379      >     0.476407,0.462092,0.447135,0.431648,0.416302,0.400932,
1380      >     0.385885,0.370871,0.356128,0.341562,0.327354,0.312933,
1381      >     0.298890,0.284864,0.270948,0.257177,0.243326,0.229321,
1382      >     0.214967,0.199986,0.183994,0.166574,0.147078,0.125313,
1383      >     0.101089,0.075065,0.048640,0.024855,0.007368,0.000073,
1384      >     0.509409,0.510722,0.511674,0.511738,0.511286,0.510766,
1385      >     0.509198,0.507564,0.505152,0.502704,0.499650,0.496452,
1386      >     0.492756,0.488747,0.483986,0.478817,0.472835,0.466164,
1387      >     0.458412,0.449362,0.438892,0.426796,0.412430,0.395391,
1388      >     0.375497,0.352110,0.324651,0.292275,0.254919,0.212461,
1389      >     0.165886,0.117521,0.071431,0.033526,0.009193,0.000091,
1390      >     0.570888,0.568118,0.564921,0.561074,0.555970,0.550478,
1391      >     0.544438,0.537826,0.530394,0.523005,0.514685,0.506342,
1392      >     0.497778,0.489102,0.479852,0.470348,0.460202,0.449558,
1393      >     0.438121,0.426046,0.412891,0.398444,0.382493,0.364765,
1394      >     0.344689,0.321925,0.295886,0.265915,0.231610,0.192817,
1395      >     0.150406,0.106352,0.064426,0.030028,0.008099,0.000077,
1396      >     0.605452,0.599834,0.593296,0.586627,0.578039,0.569200,
1397      >     0.559754,0.549624,0.538672,0.527535,0.516054,0.504396,
1398      >     0.492657,0.481104,0.468974,0.456682,0.443949,0.431100,
1399      >     0.417683,0.403750,0.389161,0.373852,0.357151,0.339250,
1400      >     0.319633,0.297684,0.273111,0.245042,0.213167,0.177272,
1401      >     0.138115,0.097498,0.058870,0.027271,0.007249,0.000067,
1402      >     0.621534,0.614157,0.605835,0.596272,0.585946,0.574335,
1403      >     0.562822,0.550233,0.536841,0.523598,0.509920,0.496293,
1404      >     0.482476,0.468781,0.455150,0.441275,0.427053,0.412946,
1405      >     0.398426,0.383702,0.368436,0.352438,0.335827,0.317985,
1406      >     0.298878,0.277758,0.254411,0.228040,0.198146,0.164659,
1407      >     0.128112,0.090261,0.054331,0.025017,0.006563,0.000059,
1408      >     0.625646,0.616575,0.607120,0.596395,0.584780,0.572447,
1409      >     0.558731,0.544391,0.529948,0.515041,0.499832,0.485344,
1410      >     0.470145,0.455555,0.440783,0.425851,0.410922,0.396061,
1411      >     0.381038,0.365825,0.350234,0.334195,0.317648,0.300202,
1412      >     0.281570,0.261329,0.238985,0.213898,0.185675,0.154079,
1413      >     0.119726,0.084159,0.050488,0.023108,0.005981,0.000052,
1414      >     0.621816,0.611958,0.601503,0.590515,0.577713,0.564667,
1415      >     0.550311,0.535181,0.519529,0.503933,0.488235,0.472641,
1416      >     0.456837,0.441572,0.426210,0.410935,0.395680,0.380387,
1417      >     0.365098,0.349710,0.334114,0.318243,0.301798,0.284744,
1418      >     0.266642,0.247149,0.225778,0.201848,0.175065,0.145114,
1419      >     0.112575,0.078976,0.047163,0.021485,0.005492,0.000047,
1420      >     0.612456,0.602615,0.592112,0.580629,0.567547,0.554000,
1421      >     0.539266,0.523505,0.507494,0.491317,0.475170,0.459335,
1422      >     0.443356,0.427588,0.412089,0.396613,0.381018,0.365836,
1423      >     0.350418,0.335168,0.319712,0.303995,0.287942,0.271269,
1424      >     0.253737,0.234846,0.214374,0.191514,0.165911,0.137420,
1425      >     0.106436,0.074522,0.044401,0.020096,0.005079,0.000042,
1426      >     0.599906,0.590658,0.579871,0.568198,0.555552,0.541752,
1427      >     0.526573,0.511088,0.494690,0.478321,0.462067,0.446064,
1428      >     0.430078,0.414174,0.398700,0.383180,0.367731,0.352392,
1429      >     0.337125,0.322012,0.306729,0.291406,0.275582,0.259369,
1430      >     0.242384,0.224158,0.204420,0.182500,0.157984,0.130669,
1431      >     0.101087,0.070603,0.041960,0.018893,0.004724,0.000039,
1432      >     0.586385,0.576832,0.566175,0.554843,0.542371,0.528570,
1433      >     0.513481,0.497995,0.481584,0.465454,0.449166,0.432955,
1434      >     0.417165,0.401569,0.386012,0.370517,0.355274,0.340119,
1435      >     0.324992,0.310107,0.295066,0.280058,0.264679,0.248876,
1436      >     0.232336,0.214712,0.195681,0.174542,0.150960,0.124776,
1437      >     0.096383,0.067222,0.039821,0.017842,0.004416,0.000035,
1438      >     0.571562,0.562230,0.553483,0.541106,0.528820,0.515154,
1439      >     0.500724,0.484959,0.468895,0.452796,0.436478,0.420645,
1440      >     0.404877,0.389321,0.373971,0.358762,0.343687,0.328774,
1441      >     0.313968,0.299216,0.284619,0.269828,0.254780,0.239398,
1442      >     0.223345,0.206321,0.187886,0.167537,0.144752,0.119523,
1443      >     0.092229,0.064212,0.037929,0.016920,0.004149,0.000033,
1444      >     0.556430,0.547716,0.538118,0.527507,0.515330,0.502167,
1445      >     0.487523,0.472378,0.456637,0.440513,0.424536,0.408900,
1446      >     0.393308,0.378059,0.362874,0.347899,0.333001,0.318366,
1447      >     0.303866,0.289409,0.275001,0.260570,0.245889,0.230899,
1448      >     0.215327,0.198784,0.180938,0.161191,0.139240,0.114867,
1449      >     0.088515,0.061527,0.036242,0.016094,0.003912,0.000030,
1450      >     0.541088,0.533114,0.524472,0.513843,0.502197,0.489243,
1451      >     0.475392,0.460194,0.444773,0.428968,0.413337,0.397800,
1452      >     0.382488,0.367387,0.352568,0.337799,0.323189,0.308771,
1453      >     0.294545,0.280383,0.266247,0.252107,0.237863,0.223231,
1454      >     0.208050,0.191998,0.174706,0.155516,0.134237,0.110630,
1455      >     0.085174,0.059101,0.034727,0.015356,0.003700,0.000028,
1456      >     0.526746,0.519306,0.510794,0.500680,0.489449,0.477026,
1457      >     0.463177,0.448622,0.433350,0.417867,0.402498,0.387171,
1458      >     0.372235,0.357466,0.342866,0.328463,0.314082,0.299986,
1459      >     0.285991,0.272106,0.258305,0.244422,0.230511,0.216246,
1460      >     0.201466,0.185822,0.169011,0.150402,0.129717,0.106819,
1461      >     0.082146,0.056897,0.033348,0.014690,0.003511,0.000026,
1462      >     0.508294,0.509537,0.510232,0.510364,0.510028,0.509066,
1463      >     0.507574,0.506096,0.503764,0.501269,0.498261,0.494940,
1464      >     0.491357,0.487139,0.482585,0.477518,0.471487,0.464685,
1465      >     0.456954,0.448138,0.437547,0.425359,0.410953,0.393893,
1466      >     0.373877,0.350282,0.322653,0.290145,0.252556,0.209889,
1467      >     0.163221,0.114885,0.069163,0.031966,0.008546,0.000081,
1468      >     0.569358,0.566803,0.562983,0.559001,0.554723,0.548955,
1469      >     0.542742,0.536150,0.528626,0.520981,0.513074,0.504898,
1470      >     0.496364,0.487653,0.478428,0.468815,0.458593,0.447996,
1471      >     0.436563,0.424615,0.411428,0.396895,0.381014,0.363131,
1472      >     0.342953,0.320173,0.294051,0.263852,0.229272,0.190367,
1473      >     0.147875,0.103902,0.062322,0.028615,0.007526,0.000069,
1474      >     0.603670,0.598438,0.591834,0.584727,0.576222,0.567680,
1475      >     0.558196,0.547591,0.536787,0.525654,0.514148,0.502844,
1476      >     0.491102,0.479283,0.467390,0.455046,0.442422,0.429358,
1477      >     0.416062,0.402212,0.387705,0.372222,0.355668,0.337678,
1478      >     0.318005,0.295922,0.271145,0.242997,0.210936,0.174968,
1479      >     0.135724,0.095201,0.057278,0.025971,0.006733,0.000060,
1480      >     0.619938,0.612660,0.603237,0.594399,0.584220,0.572945,
1481      >     0.560906,0.547974,0.535000,0.521640,0.508074,0.494503,
1482      >     0.480788,0.467052,0.453309,0.439430,0.425326,0.411154,
1483      >     0.396834,0.381957,0.366731,0.350903,0.334155,0.316435,
1484      >     0.297197,0.275976,0.252539,0.225982,0.196020,0.162415,
1485      >     0.125838,0.088092,0.052512,0.023816,0.006094,0.000053,
1486      >     0.623752,0.615162,0.604853,0.594199,0.582617,0.570319,
1487      >     0.556672,0.542512,0.527993,0.513108,0.498185,0.483234,
1488      >     0.468378,0.453629,0.438899,0.424139,0.409158,0.394331,
1489      >     0.379329,0.364110,0.348592,0.332662,0.316045,0.298521,
1490      >     0.279860,0.259548,0.237136,0.211942,0.183600,0.151947,
1491      >     0.117534,0.082104,0.048767,0.021982,0.005551,0.000047,
1492      >     0.619727,0.610281,0.600069,0.588100,0.575830,0.562394,
1493      >     0.548038,0.533038,0.517667,0.501730,0.486231,0.470564,
1494      >     0.454906,0.439606,0.424254,0.409150,0.393803,0.378587,
1495      >     0.363406,0.347957,0.332413,0.316594,0.300201,0.283135,
1496      >     0.264987,0.245411,0.223977,0.199987,0.173072,0.143020,
1497      >     0.110486,0.077018,0.045586,0.020431,0.005096,0.000042,
1498      >     0.610495,0.600588,0.589924,0.578181,0.565270,0.551795,
1499      >     0.537196,0.521624,0.505437,0.489150,0.472983,0.457321,
1500      >     0.441422,0.425796,0.410213,0.394710,0.379352,0.363971,
1501      >     0.348802,0.333422,0.317927,0.302300,0.286286,0.269619,
1502      >     0.252021,0.233168,0.212605,0.189631,0.163957,0.135368,
1503      >     0.104416,0.072632,0.042854,0.019101,0.004711,0.000038,
1504      >     0.597856,0.588522,0.577800,0.566289,0.553520,0.539436,
1505      >     0.524435,0.508812,0.492453,0.476231,0.459950,0.443865,
1506      >     0.427942,0.412263,0.396763,0.381179,0.365851,0.350734,
1507      >     0.335459,0.320164,0.305048,0.289712,0.274016,0.257736,
1508      >     0.240735,0.222511,0.202688,0.180671,0.156046,0.128701,
1509      >     0.099140,0.068812,0.040482,0.017948,0.004381,0.000035,
1510      >     0.584084,0.574752,0.564446,0.552705,0.540225,0.526223,
1511      >     0.511434,0.495863,0.479694,0.463424,0.447226,0.431058,
1512      >     0.415167,0.399560,0.384026,0.368648,0.353514,0.338375,
1513      >     0.323406,0.308441,0.293400,0.278400,0.263020,0.247257,
1514      >     0.230713,0.213051,0.193943,0.172756,0.149112,0.122832,
1515      >     0.094502,0.065475,0.038404,0.016946,0.004094,0.000032,
1516      >     0.569517,0.560385,0.550236,0.539113,0.526421,0.513222,
1517      >     0.498378,0.482790,0.466696,0.450744,0.434441,0.418603,
1518      >     0.402901,0.387328,0.372129,0.356950,0.341855,0.326972,
1519      >     0.312231,0.297539,0.282901,0.268148,0.253147,0.237801,
1520      >     0.221735,0.204645,0.186195,0.165769,0.142932,0.117672,
1521      >     0.090399,0.062530,0.036573,0.016065,0.003846,0.000029,
1522      >     0.554361,0.545724,0.536138,0.525127,0.513348,0.500011,
1523      >     0.485538,0.470383,0.454424,0.438520,0.422677,0.406888,
1524      >     0.391373,0.376154,0.360890,0.345988,0.331296,0.316666,
1525      >     0.302083,0.287727,0.273311,0.258904,0.244293,0.229265,
1526      >     0.213676,0.197142,0.179286,0.159481,0.137449,0.113024,
1527      >     0.086731,0.059884,0.034929,0.015279,0.003625,0.000027,
1528      >     0.539566,0.531375,0.522021,0.511690,0.499982,0.487119,
1529      >     0.474212,0.458150,0.442576,0.426822,0.411177,0.395600,
1530      >     0.380520,0.365481,0.350710,0.336033,0.321423,0.307074,
1531      >     0.292885,0.278746,0.264660,0.250555,0.236263,0.221649,
1532      >     0.206449,0.190350,0.173011,0.154255,0.132811,0.108865,
1533      >     0.083444,0.057506,0.033461,0.014573,0.003428,0.000025,
1534      >     0.524865,0.517319,0.508721,0.498692,0.487348,0.474698,
1535      >     0.461245,0.446504,0.431299,0.416072,0.400507,0.385278,
1536      >     0.370252,0.355605,0.341042,0.326661,0.312354,0.298301,
1537      >     0.284289,0.270473,0.256711,0.242903,0.228985,0.214684,
1538      >     0.199885,0.184233,0.167367,0.148740,0.128003,0.105093,
1539      >     0.080442,0.055358,0.032128,0.013936,0.003252,0.000024/
1540
1541 C     IF (X.GE.XMAXX.OR.Q2.LT.Q2MIN.OR.Q2.GT.Q2MAX) THEN
1542 C        WRITE(*,*) 'X or Q2 out of range!!'
1543 C        WRITE(*,*) 'Valid range is:'
1544 C        WRITE(*,*) '1E-05 < X  < 0.95 and'
1545 C        WRITE(*,*) '4     < Q2 < 520 GeV**2'
1546 C        RETURN
1547 C     END IF
1548
1549       IF (X.GE.XMAX.OR.B.GT.BMAX) THEN
1550          RES(1) = 1.0
1551          RES(2) = 1.0
1552          RETURN
1553       END IF
1554       IF(Q2.LT.Q2MIN) THEN
1555          Q2 = Q2MIN
1556       ELSE IF(Q2.GT.Q2MAX) THEN
1557          Q2 = Q2MAX
1558       ENDIF
1559
1560       CALL GGINTER(INUCL,X,Q2,SHAD)
1561
1562       DO I=1,31
1563          TATMP(I) = TA(I,INUCL)
1564       ENDDO
1565       CALL SPLINE(TMAX,IMPAR,TATMP,C,D,E)
1566       TAF = SEVAL(TMAX,B,IMPAR,TATMP,C,D,E)
1567
1568       DO IK=1,2
1569          RES(IK) = 1.0/(1.0 + (ANUCL(INUCL)-1.0)*SHAD(IK)*TAF)
1570       ENDDO
1571
1572       RETURN
1573       END
1574
1575       SUBROUTINE GGINTER(INUCL,X,Q2,SHAD)
1576       IMPLICIT DOUBLE PRECISION(A-H,L-Z)
1577       IMPLICIT INTEGER(I-K)
1578
1579       DIMENSION SHAD(2)
1580       DIMENSION G(36,13,4),LQ(36,13,4)
1581       DIMENSION XB(36),Q2V(13)
1582       DIMENSION GQTMP(13),LQQTMP(13)
1583       DIMENSION GXTMP(36),LQXTMP(36)
1584       DIMENSION C(100),D(100),E(100)
1585
1586       PARAMETER(INMAX=13,IMMAX=36)
1587
1588       COMMON/GG07/ XB,Q2V,G,LQ
1589
1590       DO K=1,36
1591          DO J=1,13
1592             GQTMP(J)  = G(K,J,INUCL)
1593             LQQTMP(J) = LQ(K,J,INUCL)
1594          ENDDO
1595          CALL SPLINE(INMAX,Q2V,GQTMP,C,D,E)
1596          GXTMP(K)  = SEVAL(INMAX,Q2,Q2V,GQTMP,C,D,E)
1597          CALL SPLINE(INMAX,Q2V,LQQTMP,C,D,E)
1598          LQXTMP(K)  = SEVAL(INMAX,Q2,Q2V,LQQTMP,C,D,E)
1599       ENDDO
1600
1601       CALL SPLINE(IMMAX,XB,GXTMP,C,D,E)
1602       SHAD(1) = SEVAL(IMMAX,X,XB,GXTMP,C,D,E)
1603       CALL SPLINE(IMMAX,XB,LQXTMP,C,D,E)
1604       SHAD(2) = SEVAL(IMMAX,X,XB,LQXTMP,C,D,E)
1605
1606       RETURN
1607       END
1608
1609 C ---------------------------------------------------------------------
1610       SUBROUTINE SPLINE(N,X,Y,B,C,D)
1611 C ---------------------------------------------------------------------
1612 c***************************************************************************
1613 C     CALCULATE THE COEFFICIENTS B,C,D IN A CUBIC SPLINE INTERPOLATION.
1614 C     INTERPOLATION SUBROUTINES ARE TAKEN FROM
1615 C     G.E. FORSYTHE, M.A. MALCOLM AND C.B. MOLER,
1616 C     COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS (PRENTICE-HALL, 1977).
1617       IMPLICIT double precision (A-H,O-Z)
1618       DIMENSION X(100),Y(100),B(100),C(100),D(100)
1619       NM1=N-1
1620       IF(N.LT.2) RETURN
1621       IF(N.LT.3) GO TO 250
1622       D(1)=X(2)-X(1)
1623       C(2)=(Y(2)-Y(1))/D(1)
1624       DO 210 I=2,NM1
1625         D(I)=X(I+1)-X(I)
1626         B(I)=2.0D0*(D(I-1)+D(I))
1627         C(I+1)=(Y(I+1)-Y(I))/D(I)
1628         C(I)=C(I+1)-C(I)
1629  210  CONTINUE
1630       B(1)=-D(1)
1631       B(N)=-D(N-1)
1632       C(1)=0.0D0
1633       C(N)=0.0D0
1634       IF(N.EQ.3) GO TO 215
1635       C(1)=C(3)/(X(4)-X(2))-C(2)/(X(3)-X(1))
1636       C(N)=C(N-1)/(X(N)-X(N-2))-C(N-2)/(X(N-1)-X(N-3))
1637       C(1)=C(1)*D(1)**2.0D0/(X(4)-X(1))
1638       C(N)=-C(N)*D(N-1)**2.0D0/(X(N)-X(N-3))
1639  215  CONTINUE
1640       DO 220 I=2,N
1641         T=D(I-1)/B(I-1)
1642         B(I)=B(I)-T*D(I-1)
1643         C(I)=C(I)-T*C(I-1)
1644  220  CONTINUE
1645       C(N)=C(N)/B(N)
1646       DO 230 IB=1,NM1
1647         I=N-IB
1648         C(I)=(C(I)-D(I)*C(I+1))/B(I)
1649  230  CONTINUE
1650       B(N)=(Y(N)-Y(NM1))/D(NM1)+D(NM1)*(C(NM1)+2.0D0*C(N))
1651       DO 240 I=1,NM1
1652         B(I)=(Y(I+1)-Y(I))/D(I)-D(I)*(C(I+1)+2.0D0*C(I))
1653         D(I)=(C(I+1)-C(I))/D(I)
1654         C(I)=3.0D0*C(I)
1655  240  CONTINUE
1656       C(N)=3.0D0*C(N)
1657       D(N)=D(N-1)
1658       RETURN
1659  250  CONTINUE
1660       B(1)=(Y(2)-Y(1))/(X(2)-X(1))
1661       C(1)=0.0D0
1662       D(1)=0.0D0
1663       B(2)=B(1)
1664       C(2)=0.0D0
1665       D(2)=0.0D0
1666       RETURN
1667       END
1668 c
1669 c***************************************************************************
1670 C ---------------------------------------------------------------------
1671       double precision FUNCTION SEVAL(N,XX,X,Y,B,C,D)
1672 C ---------------------------------------------------------------------
1673 c***************************************************************************
1674 C CALCULATE THE DISTRIBUTION AT XX BY CUBIC SPLINE INTERPOLATION.
1675       implicit double precision(A-H,O-Z)
1676       DIMENSION X(100),Y(100),B(100),C(100),D(100)
1677       DATA I/1/
1678       IF(I.GE.N) I=1
1679       IF(XX.LT.X(I)) GO TO 310
1680       IF(XX.LE.X(I+1)) GO TO 330
1681  310  CONTINUE
1682       I=1
1683       J=N+1
1684  320  CONTINUE
1685       K=(I+J)/2
1686       IF(XX.LT.X(K)) J=K
1687       IF(XX.GE.X(K)) I=K
1688       IF(J.GT.I+1) GO TO 320
1689  330  CONTINUE
1690       DX=XX-X(I)
1691       SEVAL=Y(I)+DX*(B(I)+DX*(C(I)+DX*D(I)))
1692       RETURN
1693       END
1694 c
1695 c***************************************************************************
1696
1697
1698 * function to generate gauss distribution
1699       double precision function gaus(x0,sig)
1700       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1701  41   u1=pyr(0)
1702       u2=pyr(0)
1703       v1=2.d0*u1-1.d0
1704       v2=2.d0*u2-1.d0
1705       s=v1**2+v2**2
1706       if(s.gt.1) go to 41
1707       gaus=v1*dsqrt(-2.d0*dlog(s)/s)*sig+x0
1708       return
1709       end