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