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