*---------------------------------------------------------------------- * * Filename : HYDJET1_8.F * * Author : Igor Lokhtin (Igor.Lokhtin@cern.ch) * Version : HYDJET1_8.f * Last revision : 17-JUL-2011 * *====================================================================== * * Description : Event generator for simulation of jet production, jet * quenching and flow effects in ultrarelativistic heavy * ion AA collisions * * Method: I.P. Lokhtin and A.M. Snigirev, Eur. Phys. J C 45 (2006) 211 * *====================================================================== ********************************* HYINIT *************************** SUBROUTINE HYINIT(energy,A,ifb1,bmin,bmax,bfix1,nh1) * PYTHIA inizialization, calculation of total and hard cross sections and * # of participants and binary sub-collisions at "reference point" (Pb,b=0), * tabulation of nuclear thickness function and nuclear overlap function, * test of input parametes IMPLICIT DOUBLE PRECISION(A-H, O-Z) double precision numpar,npar0,nbco0 external numpar,rhoaa,hfunc3 common /pyint7/ sigt(0:6,0:6,0:5) common /pypars/ mstp(200),parp(200),msti(200),pari(200) common /pysubs/ msel,mselpd,msub(500),kfin(2,-40:40),ckin(200) common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr common /hypyin/ ene,rzta,rnta,bfix,ifb,nh common /hyflow/ ytfl,ylfl,Tf,fpart common /hygeom/ BC common /hythic/ BAB(110),TAB(110),TAAB(110) common /hynup1/ bp,x save /pyint7/,/pypars/,/pysubs/,/hyjpar/,/hyipar/,/hypyin/, > /hyflow/,/hygeom/,/hythic/ * start HYDJET initialization init=1 ipr=1 * set beam paramters ene=energy ! c.m.s. energy per nucleon AW=A ! atomic weight RA=1.15d0*AW**0.333333d0 ! nuclear radius rzta=1.d0/(1.98d0+0.015d0*A**0.666667d0) ! fraction of protons in nucleus rnta=1.d0-rzta ! fraction of neutrons in nucleus ifb=ifb1 ! centrality flag bfix=bfix1 ! fixed impact parameter nh=nh1 ! mean soft mult. in central PbPb ptmin=ckin(3) ! minimum pt of hard scattering * printout of HYDJET versioning write(6,*) '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >%%%%%%%%%%%%%%%%%%%%%%%%' write(6,*) '%% This is HYDJET version 1.8 > %%' write(6,*) '%% Corresponding author: Igor Lokhtin (Igor.Lo >khtin@cern.ch) %%' write(6,*) '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >%%%%%%%%%%%%%%%%%%%%%%%%' * Pythia inizialization for pp collisions call pyinit('cms','p','p',ene) * calculation of hard scattering cross section if(nhsel.ne.0) then mstp(111)=0 * no printout of Pythia initialization information hereinafter mstp(122)=0 * Pythia test pp event run do i=1,1000 call pyevnt end do * hard scattering pp cross section sjpp=pari(1) * Pythia inizialization for pn collisions call pyinit('cms','p','n',ene) * Pythia test pn event run do i=1,1000 call pyevnt end do * hard scattering pn cross section sjpn=pari(1) * Pythia inizialization for nn collisions call pyinit('cms','n','n',ene) * Pythia test nn event run do i=1,1000 call pyevnt end do * hard scattering nn cross section sjnn=pari(1) sigjet=rzta*rzta*sjpp+rnta*rnta*sjnn+2.d0*rzta*rnta*sjpn end if * total inelastic cross section if(sigin.lt.10.d0.or.sigin.gt.200.d0) > sigin=sigt(0,0,0)-sigt(0,0,1) * # of nucelons-participants and NN sub-collisions at "reference point" (Pb,b=0) Apb=207.d0 Rpb=1.15d0*Apb**0.333333d0 EPS=0.005d0 Z2=4.d0*Rpb Z1=-1.d0*Z2 H=0.01d0*(Z2-Z1) do ib=1,110 BC=3.d0*Rpb*(ib-1)/109.d0 CALL SIMPA(Z1,Z2,H,EPS,1.d-8,rhoaa,Z,RES,AIH,AIABS) BAB(ib)=BC TAB(ib)=Apb*RES end do Z1=0.d0 Z2=6.28318d0 H=0.01d0*(Z2-Z1) bp=0.d0 CALL SIMPA(Z1,Z2,H,EPS,1.d-8,HFUNC3,X,TAAPB0,AIH,AIABS) npar0=numpar(0.d0) ! Npart(Pb,b=0) nbco0=0.1d0*sigin*TAAPB0 ! Nsub(Pb,b=0) init=0 * creation of arrays for tabulation of beam/target nuclear thickness function Z2=4.d0*RA Z1=-1.d0*Z2 H=0.01d0*(Z2-Z1) do ib=1,110 BC=3.d0*RA*(ib-1)/109.d0 CALL SIMPA(Z1,Z2,H,EPS,1.d-8,rhoaa,Z,RES,AIH,AIABS) BAB(ib)=BC TAB(ib)=AW*RES end do * creation of arrays for tabulation of nuclear overlap function Z1=0.d0 Z2=6.28318d0 H=0.01d0*(Z2-Z1) do ib=1,110 bp=BAB(ib) CALL SIMPA(Z1,Z2,H,EPS,1.d-8,HFUNC3,X,RES,AIH,AIABS) TAAB(ib)=RES end do * test of centrality selection if(ifb.eq.0) then if(bfix.lt.0.d0) then write(6,*) 'Impact parameter less than zero!' bfix=0.d0 end if if (bfix.gt.3.d0) then write(6,*) 'Impact parameter larger than three nuclear radius!' bfix=3.d0 end if else if(bmin.lt.0.d0) then write(6,*) 'Impact parameter less than zero!' bmin=0.d0 end if if(bmax.gt.3.d0) then write(6,*) 'Impact parameter larger than three nuclear radius!' bmax=3.d0 end if bminh=bmin bmaxh=bmax end if * test of flow parameter selection if (Tf.lt.0.08d0.or.Tf.gt.0.2d0) Tf=0.1d0 ! freeze-out temperature if (ylfl.lt.0.01d0.or.ylfl.gt.7.d0) ylfl=4.d0 ! longitudinal flow rapidity if (ytfl.lt.0.01d0.or.ytfl.gt.3.d0) ytfl=1.5d0 ! transverse flow rapidity if (fpart.le.0.d0.or.fpart.gt.1.d0) fpart=1.d0 ! fraction of soft multiplicity ! proport. to # of participants * test of 'nhsel' selection if(nhsel.ne.1.and.nhsel.ne.2.and.nhsel.ne.3.and.nhsel.ne.4) > nhsel=0 return end ****************************** END HYINIT *************************** ********************************* HYEVNT **************************** SUBROUTINE HYEVNT * generation of single HYDJET event (soft+hard parts) at given parameters IMPLICIT DOUBLE PRECISION(A-H, O-Z) double precision numpar,npar0,nbco0,npart,nbcol real p,v INTEGER PYK,PYCHGE,PYCOMP external hsin,gauss,hftaa,numpar,hyhard,hipsear,pyr,pymass external ludata common /pyjets/ nhi,npadhi,khi(4000,5),phi(4000,5),vhi(4000,5) COMMON /PYDAT3/ MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5) common /lujets/ n,k(150000,5),p(150000,5),v(150000,5) common /hyjets/ nhj,nhp,khj(150000,5),phj(150000,5),vhj(150000,5) common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr common /hypyin/ ene,rzta,rnta,bfix,ifb,nh common /hyfpar/ bgen,nbcol,npart,npyt,nhyd common /hyflow/ ytfl,ylfl,Tf,fpart common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet save /pyjets/,/pydat3/,/lujets/,/hyjets/,/hyipar/,/hyfpar/, + /hyflow/,/hyjpar/,/hypyin/ * reset lujets and hyjets arrays before event generation n=0 nhj=0 do ncl=1,150000 do j=1,5 p(ncl,j)=0. phj(ncl,j)=0.d0 v(ncl,j)=0. vhj(ncl,j)=0.d0 k(ncl,j)=0 khj(ncl,j)=0 enddo end do * pi=3.14159d0 * generate impact parameter of A-A collision if(ifb.eq.0) then b1=bfix*RA bgen=bfix else call hipsear(fmax1,xmin1) fmax=fmax1 xmin=xmin1 3 bb1=xmin*pyr(0)+bminh*RA ff1=fmax*pyr(0) fb=hsin(bb1) if(ff1.gt.fb) goto 3 b1=bb1 bgen=bb1/RA end if * calculate # of nucelons-participants and binary NN sub-collisions npart=numpar(b1) ! Npart(b) nbcol=0.1d0*sigin*hftaa(b1) ! Nsub(b) * generate soft multiplicity in event, np, * fpart - fraction of soft multiplicity proportional to # of participants, * fbcol=1-fpart - fraction of multiplicity proprtional to # of NN subcollisions fbcol=1.d0-fpart rnp=nh*(fpart*npart+fbcol*nbcol)/(fpart*npar0+fbcol*nbco0) np=int(rnp) sign=dsqrt(rnp) 1 if(nhsel.lt.3.and.np.gt.0) then np=max(0,int(gauss(rnp,sign))) else np=0 end if if(np.gt.150000) then write(6,*) 'Warning, soft multiplicity too large!' goto 1 end if * generate hard parton-parton scatterings (Q>ptmin) 'njet' times njet=0 if(nhsel.ne.0) then pjet=sigjet/sigin do i=1,int(nbcol) if(pyr(0).lt.pjet) njet=njet+1 end do call hyhard end if * reset pyjets array before soft component generation nhi=0 do ncl=1,4000 do j=1,5 phi(ncl,j)=0.d0 vhi(ncl,j)=0.d0 khi(ncl,j)=0 enddo end do npyt=nhj nhyd=np if(nhyd.eq.0) goto 4 if(bgen.gt.1.999d0) b1=1.999d0*RA * set of parameters for azimuthal anisotropy in hadron emission region anc=2.8d0 ank=0.25d0 ane=0.5d0*ank*b1/RA anb=max(0.000000000001d0,anc*(1.-ane*ane)*ane) anl=max(0.000000000001d0,0.5d0*(dsqrt(1.d0+4.d0*anb* > (ane+anb))-1.d0)/anb) anr=RA*dsqrt(abs(1.d0-ane)) * generate sort of hadrons in jetset/pythia format 5 yy=9999.d0*pyr(0) if(yy.lt.1376.d0) then kf=211 elseif(yy.lt.2752.d0) then kf=-211 elseif(yy.lt.4142.d0) then kf=111 elseif(yy.lt.4590.d0) then kf=321 elseif(yy.lt.5038.d0) then kf=-321 elseif(yy.lt.5479.5d0) then kf=311 elseif(yy.lt.5921.d0) then kf=-311 elseif(yy.lt.6063.5d0) then kf=2212 elseif(yy.lt.6206.d0) then kf=-2212 elseif(yy.lt.6347.7d0) then kf=2112 elseif(yy.lt.6489.4d0) then kf=-2112 elseif(yy.lt.6935.4d0) then kf=213 elseif(yy.lt.7381.4d0) then kf=-213 elseif(yy.lt.7827.4d0) then kf=113 elseif(yy.lt.8248.8d0) then kf=223 elseif(yy.lt.8614.5d0) then kf=221 elseif(yy.lt.8680.d0) then kf=331 elseif(yy.lt.8829.d0) then kf=333 elseif(yy.lt.9092.d0) then kf=323 elseif(yy.lt.9355.d0) then kf=-323 elseif(yy.lt.9613.d0) then kf=313 elseif(yy.lt.9871.d0) then kf=-313 elseif(yy.lt.9935.d0) then kf=3122 else kf=-3122 end if nhj=nhj+1 khj(nhj,1)=1 khj(nhj,2)=kf do j=3,5 khj(nhj,j)=0 end do do j=1,5 vhj(nhj,j)=0.d0 enddo kfa=iabs(kf) phj(nhj,5)=pymass(kfa) * generate uniform distribution in system of a fluid element 2 ep0=-1.d0*Tf*(dlog(max(1.d-10,pyr(0)))+dlog(max(1.d-10,pyr(0))) > +dlog(max(1.d-10,pyr(0)))) if(ep0.le.phj(nhj,5)) go to 2 pp0=dsqrt(abs(ep0**2-abs(phj(nhj,5)**2))) probt=pp0/ep0 if(pyr(0).gt.probt) go to 2 ctp0=2.d0*pyr(0)-1.d0 stp0=dsqrt(abs(1.d0-ctp0**2)) php0=2.d0*pi*pyr(0) * generate coordinates of a fluid element c etaf=ylfl*(2.d0*pyr(0)-1.d0) ! flat initial eta-spectrum etaf=gauss(0.d0,ylfl) ! gaussian initial eta-spectrum phif=2.d0*pi*pyr(0) RF=anr*dsqrt(abs(1.d0-ane*ane))/ > (dsqrt(abs(1.d0+ane*dcos(2.d0*phif)))) rrf=RF*RF*pyr(0) * generate four-velocity of a fluid element vradf=dsinh(ytfl*dsqrt(rrf)/anr) rn1=dcosh(etaf) rn2=dsinh(etaf) uxf=vradf*dsqrt(abs(1.d0+anl))*dcos(phif) uyf=vradf*dsqrt(abs(1.d0-anl))*dsin(phif) urf=dsqrt(uxf*uxf+uyf*uyf) rn0=dsqrt(abs(1.d0+urf*urf)) utf=rn1*rn0 uzf=rn2*rn0 * Von Neumann rejection/acceptance procedure for Cooper-Frye prescription rfac=(rn1+rn0)/(1.d0+utf) rnx=-1.d0*uxf*rfac rny=-1.d0*uyf*rfac rnz=rn2-uzf*rfac pnx=pp0*stp0*dcos(php0) pny=pp0*stp0*dsin(php0) pnz=pp0*ctp0 pn0=dsqrt(pnx**2+pny**2+pnz**2+phj(nhj,5)**2) wesn=rn0-(rnx*pnx+rny*pny+rnz*pnz)/pn0 wksi=2.d0*dcosh(ytfl)*pyr(0) if(wksi.gt.wesn) goto 2 * boost in the laboratory system uipi=pp0*(urf*stp0*dcos(phif-php0)+uzf*ctp0) bfac=uipi/(utf+1.d0)+ep0 phj(nhj,1)=pnx+urf*dcos(phif)*bfac phj(nhj,2)=pny+urf*dsin(phif)*bfac phj(nhj,3)=pnz+uzf*bfac phj(nhj,4)=dsqrt(phj(nhj,1)**2+phj(nhj,2)**2+phj(nhj,3)**2 > +phj(nhj,5)**2) * treat hadron decay (if it switched on) with pyexec routine if(mdcy(pycomp(kf),1).eq.1) then nhi=1 do j=1,5 phi(1,j)=phj(nhj,j) khi(1,j)=khj(nhj,j) vhi(1,j)=vhj(nhj,j) end do call pyexec do i=1,nhi do j=1,5 icor=0 if(j.gt.2.and.khi(i,j).gt.0) icor=nhj-1 phj(nhj+i-1,j)=phi(i,j) khj(nhj+i-1,j)=khi(i,j)+icor vhj(nhj+i-1,j)=vhi(i,j) end do end do nhj=nhj+nhi-1 np=np+nhi-1 end if if(nhj.lt.npyt+np) goto 5 nhyd=np 4 continue * fill array 'lujets' n=nhyd+npyt do ih=1,nhj do jh=1,5 p(ih,jh)=phj(ih,jh) k(ih,jh)=khj(ih,jh) v(ih,jh)=vhj(ih,jh) end do end do return end ****************************** END HYEVNT ************************** ********************************* HYHARD *************************** SUBROUTINE HYHARD * generate 'njet' number of hard parton-parton scatterings with PYTHIA IMPLICIT DOUBLE PRECISION(A-H, O-Z) double precision npar0,nbco0,npart,nbcol INTEGER PYK,PYCHGE,PYCOMP CHARACTER beam*2,targ*2 external pydata external pyp,pyr,pyk,pyquen,shad1 common /pyjets/ n,npad,k(4000,5),p(4000,5),v(4000,5) common /hyjets/ nhj,nhp,khj(150000,5),phj(150000,5),vhj(150000,5) COMMON /PYDAT1/ MSTU(200),PARU(200),MSTJ(200),PARJ(200) COMMON /PYDAT2/ KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4) COMMON /PYDAT3/ MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5) COMMON /PYSUBS/ MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200) common /pypars/ mstp(200),parp(200),msti(200),pari(200) common /parimp/ b1, psib1, rb1, rb2, noquen common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr common /hyfpar/ bgen,nbcol,npart,npyt,nhyd common /hypyin/ ene,rzta,rnta,bfix,ifb,nh save /pyjets/,/pypars/,/pydat1/,/pydat2/,/pydat3/,/pysubs/, + /hyjets/,/parimp/,/hyjpar/,/hyipar/,/hyfpar/,/hypyin/ * generate 'njet' PYTHIA events and fill arrays for partons and hadrons nshad=0 noquen=0 if(nhsel.eq.1.or.nhsel.eq.3) noquen=1 if(njet.ge.1) then ifbp=0 ! fix impact parameter bfixp=RA*bgen Ap=AW ! atomic weight do ihard=1,njet mstp(111)=0 * generate type of NN sub-collision (pp, pn or nn) rand1=pyr(0) if(rand1.lt.rzta) then beam='p' else beam='n' end if rand2=pyr(0) if(rand2.lt.rzta) then targ='p' else targ='n' end if call pyinit('cms',beam,targ,ene) c mstj(41)=0 ! vacuum showering off call pyevnt ! generate hard scattering * PYQUEN: quenched jets if noquen=0 or non-quenched jets if noquen=1 and ishad=1 if(ishad.eq.1.or.nhsel.eq.2.or.nhsel.eq.4) > call pyquen(Ap,ifbp,bfixp) * treatment of "nuclear shadowing" (for Pb, Au, Pd or Ca beams only) if(ishad.eq.1) then kfh1=abs(k(3,2)) kfh2=abs(k(4,2)) xh1=pari(33) xh2=pari(34) Q2=pari(22) shad=shad1(kfh1,xh1,Q2,rb1)*shad1(kfh2,xh2,Q2,rb2) if(pyr(0).gt.shad) then nshad=nshad+1 goto 53 end if end if call pyexec ! hadronization done c call pyedit(2) ! remove partons & leave hadrons * fill array of final particles nu=nhj+n if(nu.gt.150000-np) then write(6,*) 'Warning, multiplicity too large! Cut hard part.' goto 52 end if nhj=nu do i=nhj-n+1,nhj ip=i+n-nhj do j=1,5 phj(i,j)=p(ip,j) end do do j=1,5 vhj(i,j)=v(ip,j) end do do j=1,5 khj(i,j)=k(ip,j) end do do j=3,5 kk=khj(i,j) if(kk.gt.0) then khj(i,j)=kk+nhj-n end if end do end do 53 continue end do 52 njet=ihard-1 end if njet=njet-nshad return end ****************************** END HYHARD ************************** ********************************* HIPSEAR *************************** SUBROUTINE HIPSEAR (fmax,xmin) * find maximum and 'sufficient minimum' of differential inelasic AA cross * section as a function of impact paramater (xm, fm are outputs) IMPLICIT DOUBLE PRECISION(A-H, O-Z) double precision npar0,nbco0 external hsin common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr save /hyipar/ xmin=(bmaxh-bminh)*RA fmax=0.d0 do j=1,1000 x=bminh*RA+xmin*(j-1)/999.d0 f=hsin(x) if(f.gt.fmax) then fmax=f endif end do return end ****************************** END HIPSEAR ************************** * differential inelastic AA cross section double precision function hsin(x) IMPLICIT DOUBLE PRECISION(A-H, O-Z) external hftaa common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet save /hyjpar/ br=x hsin=br*(1.d0-dexp(-0.1d0*hftaa(br)*sigin)) return end * number of nucleons-participants at impact parameter b double precision function numpar(c) IMPLICIT DOUBLE PRECISION(A-H, O-Z) external HFUNC1 common /hynup1/ bp,x EPS=0.005d0 A=0.d0 B=6.28318d0 H=0.01d0*(B-A) bp=c CALL SIMPA(A,B,H,EPS,1.d-8,HFUNC1,X,RES,AIH,AIABS) numpar=2.d0*RES return end * double precision function HFUNC1(x) IMPLICIT DOUBLE PRECISION(A-H, O-Z) double precision npar0,nbco0 external HFUNC2 common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr common /hynup1/ bp,xx save /hyipar/ if(init.eq.1) then Rl=Rpb else Rl=RA end if xx=x EPS=0.005d0 A=0.d0 B=3.d0*Rl H=0.01d0*(B-A) CALL SIMPB(A,B,H,EPS,1.d-8,HFUNC2,Y,RES,AIH,AIABS) HFUNC1=RES return end * double precision function HFUNC2(y) IMPLICIT DOUBLE PRECISION(A-H, O-Z) double precision npar0,nbco0 external hythik common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet common /hynup1/ bp,x save /hyipar/,/hyjpar/ r1=dsqrt(abs(y*y+bp*bp/4.d0+y*bp*dcos(x))) r2=dsqrt(abs(y*y+bp*bp/4.d0-y*bp*dcos(x))) s=1.d0-dexp(-0.1d0*sigin*hythik(r2)) HFUNC2=y*hythik(r1)*s return end * nuclear overlap function at impact parameter b double precision function hftaa(c) IMPLICIT DOUBLE PRECISION(A-H, O-Z) common /hythic/ BAB(110),TAB(110),TAAB(110) save /hythic/ call parinv(c,BAB,TAAB,110,RES) hftaa=RES return end * double precision function HFUNC3(x) IMPLICIT DOUBLE PRECISION(A-H, O-Z) double precision npar0,nbco0 external HFUNC4 common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr common /hynup1/ bp,xx save /hyipar/ if(init.eq.1) then Rl=Rpb else Rl=RA end if xx=x EPS=0.005d0 A=0.d0 B=3.d0*Rl H=0.01d0*(B-A) CALL SIMPB(A,B,H,EPS,1.d-8,HFUNC4,Y,RES,AIH,AIABS) HFUNC3=RES return end * double precision function HFUNC4(y) IMPLICIT DOUBLE PRECISION(A-H, O-Z) double precision npar0,nbco0 external hythik common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet common /hynup1/ bp,x save /hyipar/,/hyjpar/ r1=dsqrt(abs(y*y+bp*bp/4.d0+y*bp*dcos(x))) r2=dsqrt(abs(y*y+bp*bp/4.d0-y*bp*dcos(x))) HFUNC4=y*hythik(r1)*hythik(r2) return end * nuclear thickness function double precision function hythik(r) IMPLICIT DOUBLE PRECISION(A-H, O-Z) common /hythic/ BAB(110),TAB(110),TAAB(110) save /hythic/ call parinv(r,BAB,TAB,110,RES) hythik=RES return end * Wood-Saxon nucleon distrubution double precision function rhoaa(z) IMPLICIT DOUBLE PRECISION(A-H, O-Z) double precision npar0,nbco0 common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr common /hygeom/ BC save /hyipar/,/hygeom/ if(init.eq.1) then Rl=Rpb else Rl=RA end if pi=3.14159d0 df=0.54d0 r=sqrt(bc*bc+z*z) rho0=3.d0/(4.d0*pi*Rl**3)/(1.d0+(pi*df/Rl)**2) rhoaa=rho0/(1.d0+dexp((r-Rl)/df)) return end * function to calculate nuclear shadowing factor (for Pb, Au, Pd or Ca beams) double precision function shad1(kfh,xbjh,Q2h,r) IMPLICIT DOUBLE PRECISION(A-H, O-Z) double precision npar0,nbco0,nbcol,npart external ggshad common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr common /hyfpar/ bgen,nbcol,npart,npyt,nhyd common /hyshad/ bbmin,bbmax,inuc save /hyipar/,/hyfpar/,/hyshad/ dimension res(2) kf=kfh xbj=xbjh Q2=Q2h bb=r inuc=0 if(AW.gt.205.d0.and.AW.lt.209.d0) inuc=4 ! Pb-206, 207 or 208 if(AW.gt.196.d0.and.AW.lt.198.d0) inuc=3 ! Au-197 if(AW.gt.109.d0.and.AW.lt.111.d0) inuc=2 ! Pd-110 if(AW.gt.39.d0.and.AW.lt.41.d0) inuc=1 ! Ca-40 if(inuc.eq.0.and.ipr.eq.1) then write(6,*) > 'Warning! Shadowing is not foreseen for atomic weigth A =' > ,AW write(6,*)'****************************************************** >************************' ipr=0 end if if(inuc.ne.0) then xbj=max(5.d-5,xbj) xbj=min(0.95d0,xbj) Q2=max(4.d0,Q2) Q2=min(520.d0,Q2) bb=max(0.d0,bb) call ggshad(inuc,xbj,Q2,bb,res,ta) if(kf.eq.21) then shad1=res(1) elseif(kf.eq.1.or.kf.eq.2.or.kf.eq.3) then shad1=res(2) else shad1=1.d0 end if else shad1=1.d0 end if return end ****************************************************************************** * The part of the code which follows below, includes nuclear shadowing model * * and has been written by Konrad Tywoniuk (Oslo University, Norway) * ****************************************************************************** c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$ c$$$ c$$$ c$$$ Shadowing from Glauber-Gribov theory for c$$$ gluons and light quarks (u,d,s and their antiquarks). c$$$ All Pomeron tree diagrams have been summed (Schwimmer model). c$$$ More details about the model in c$$$ K. Tywoniuk, I.C. Arsene, L. Bravina, A. Kaidalov and E. c$$$ Zabrodin, Phys. Lett. B 657 (2007) 170 c$$$ c$$$ We use FIT B from the H1 parameterization published in c$$$ A. Aktas et al. (H1 Collaboration), Eur. Phys. J C 48 (2006) 715 c$$$ A. Aktas et al. (H1 Collaboration), Eur. Phys. J C 48 (2006) 749 c$$$ c$$$ Main routine is GGSHAD which provides the user with the c$$$ ratio of nuclear and nucleon parton distribution function c$$$ normalized by the atomic number for a given c$$$ INUCL: nucleus (1=Ca,2=Pd,3=Au,4=Pb) c$$$ X: Bjorken x c$$$ Q2: q**2, scale squared c$$$ B: impact parameter c$$$ c$$$ Then RES(1): gluon shadowing c$$$ RES(2): sea quark shadowing c$$$ TA: nuclear profile function c$$$ c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$ SUBROUTINE GGSHAD(INUCL,X,Q2,B,RES,TAF) IMPLICIT NONE DOUBLE PRECISION TA(31,4),IMPAR(31),ANUCL(4) DOUBLE PRECISION XB(36),Q2V(13) DOUBLE PRECISION XMAX,XMAXX,Q2MIN,Q2MAX,BMAX DOUBLE PRECISION G(36,13,4),LQ(36,13,4) DOUBLE PRECISION TATMP(31),SHAD(2) DOUBLE PRECISION C(100),D(100),E(100) DOUBLE PRECISION X,Q2,B DOUBLE PRECISION RES(2) DOUBLE PRECISION TAF DOUBLE PRECISION SEVAL INTEGER INUCL,IREAD,TMAX INTEGER I,J,K,IK PARAMETER(TMAX=31) PARAMETER(XMAX=0.1) PARAMETER(XMAXX=0.95) PARAMETER(Q2MIN=4.) PARAMETER(Q2MAX=520.) PARAMETER(BMAX=30.) COMMON/GG07/XB,Q2V,G,LQ DATA ANUCL/40.,110.,197.,206./ DATA IMPAR > /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., > 8.5,9.,9.5,10.,10.5,11.,11.5,12.,12.5,13.,13.5,14.,14.5,15./ DATA TA > /0.0291662,0.0288636,0.0279354,0.0263044,0.0238341,0.0203565, > 0.0158333,0.0107343,0.00617758,0.00307859,0.00139738, > 0.000603839,0.000254843,0.000106329,4.40965e-05,1.82206e-05, > 7.50949e-06,3.08881e-06,1.26837e-06,5.20082e-07,2.12978e-07, > 8.71154e-08,3.55956e-08,1.45304e-08,5.92622e-09,2.41504e-09, > 9.83429e-10,4.00184e-10,1.62741e-10,6.61407e-11,2.68657e-11, > 0.0153487,0.0152775,0.0150614,0.0146926,0.0141561,0.0134268, > 0.0124651,0.0112115,0.00959258,0.00756907,0.00527514, > 0.00313482,0.00159852,0.000732479,0.000316628,0.000133116, > 5.52508e-05,2.27904e-05,9.3692e-06,3.84349e-06,1.57423e-06, > 6.43962e-07,2.63132e-07,1.07414e-07,4.38089e-08,1.78529e-08, > 7.2699e-09,2.95832e-09,1.20304e-09,4.8894e-10,1.98603e-10, > 0.0105299,0.010498,0.0104017,0.010239,0.010006,0.00969677, > 0.00930205,0.00880781,0.00819264,0.00742455,0.00646058, > 0.00526252,0.00385835,0.00243988,0.00131252,0.000621079, > 0.000272322,0.000114993,4.77312e-05,1.96569e-05,8.06377e-06, > 3.30066e-06,1.34904e-06,5.50751e-07,2.24633e-07,9.15434e-08, > 3.72777e-08,1.51694e-08,6.16885e-09,2.50714e-09,1.01838e-09, > 0.0102179,0.0101879,0.0100975,0.00994467,0.00972606, > 0.00943628,0.0090671,0.00860605,0.00803416,0.00732289, > 0.00643266,0.00532296,0.00400025,0.00261264,0.00144993, > 0.00070098,0.000310867,0.000131949,5.48887e-05,2.26247e-05, > 9.28457e-06,3.80091e-06,1.55358e-06,6.34273e-07,2.58701e-07, > 1.05427e-07,4.29316e-08,1.74701e-08,7.10447e-09,2.88739e-09, > 1.17283e-09/ DATA Q2V > /4.000, 6.000, 9.000, 13.500, 20.250, 30.375, 45.562, > 68.344, 102.516, 153.773, 230.660, 345.990, 518.985/ DATA XB > /0.99999998E-05, 0.13000000E-04, 0.16899999E-04, > 0.21970000E-04, 0.28561000E-04, 0.37129299E-04, > 0.48268099E-04, 0.62748499E-04, 0.81573096E-04, > 0.10604500E-03, 0.13785800E-03, 0.17921600E-03, > 0.23298099E-03, 0.30287501E-03, 0.39373801E-03, > 0.51185902E-03, 0.66541700E-03, 0.86504198E-03, > 0.11245500E-02, 0.14619200E-02, 0.19005000E-02, > 0.24706500E-02, 0.32118401E-02, 0.41753901E-02, > 0.54280101E-02, 0.70564100E-02, 0.91733299E-02, > 0.11925300E-01, 0.15502900E-01, 0.20153800E-01, > 0.26200000E-01, 0.34059901E-01, 0.44277899E-01, > 0.57561301E-01, 0.74829698E-01, 0.97278602E-01/ DATA G > /2.542140,2.457280,2.372200,2.283920,2.195320,2.105540, > 2.014040,1.921480,1.827600,1.733930,1.639990,1.546170, > 1.453220,1.361040,1.268750,1.178500,1.088730,1.000740, > 0.914606,0.830842,0.749315,0.671419,0.596883,0.526173, > 0.460067,0.398489,0.341895,0.289998,0.242628,0.199405, > 0.159625,0.122521,0.087252,0.053699,0.022804,0.000414, > 1.825600,1.768440,1.712550,1.654660,1.595830,1.533980, > 1.472860,1.409910,1.346060,1.282030,1.217520,1.153050, > 1.089030,1.025930,0.962709,0.900094,0.837787,0.776733, > 0.716072,0.656570,0.598295,0.541607,0.486661,0.433846, > 0.383421,0.335585,0.290685,0.248542,0.209305,0.172614, > 0.138166,0.105552,0.074246,0.044592,0.017935,0.000291, > 1.445600,1.405920,1.364150,1.322050,1.276900,1.232020, > 1.184670,1.136620,1.086950,1.037370,0.987225,0.937488, > 0.888244,0.839061,0.790370,0.741798,0.693426,0.645511, > 0.598298,0.551516,0.505500,0.460400,0.416326,0.373488, > 0.332281,0.292614,0.254859,0.218944,0.184976,0.152759, > 0.122077,0.092747,0.064477,0.037873,0.014527,0.000216, > 1.209700,1.179300,1.146560,1.112700,1.077790,1.041400, > 1.003150,0.964598,0.923825,0.882423,0.841304,0.800381, > 0.759486,0.718947,0.678587,0.638579,0.598671,0.558978, > 0.519524,0.480617,0.442116,0.404231,0.367057,0.330599, > 0.295291,0.261100,0.228220,0.196605,0.166367,0.137331, > 0.109496,0.082677,0.056835,0.032696,0.012013,0.000165, > 1.049390,1.024810,0.998839,0.971705,0.942633,0.912223, > 0.880504,0.847236,0.812716,0.777477,0.741670,0.706578, > 0.671021,0.636183,0.601384,0.566923,0.532196,0.497956, > 0.463884,0.430069,0.396612,0.363536,0.330915,0.298900, > 0.267627,0.237155,0.207697,0.179169,0.151566,0.124961, > 0.099266,0.074436,0.050544,0.028462,0.010013,0.000127, > 0.943086,0.922388,0.900700,0.877793,0.853134,0.826808, > 0.798937,0.769284,0.738709,0.706829,0.674818,0.643041, > 0.611250,0.580168,0.548982,0.517908,0.486720,0.456021, > 0.425312,0.394899,0.364616,0.334682,0.305173,0.276042, > 0.247465,0.219584,0.192474,0.166033,0.140391,0.115481, > 0.091347,0.068001,0.045636,0.025180,0.008510,0.000101, > 0.858923,0.841883,0.823827,0.803952,0.782124,0.759447, > 0.734659,0.708141,0.680281,0.651690,0.622372,0.593451, > 0.564610,0.536045,0.507521,0.479248,0.450845,0.422731, > 0.394687,0.366774,0.339095,0.311580,0.284350,0.257527, > 0.231101,0.205195,0.179929,0.155178,0.131077,0.107555, > 0.084712,0.062629,0.041557,0.022497,0.007326,0.000082, > 0.790086,0.776082,0.760408,0.743282,0.724441,0.703812, > 0.681673,0.657893,0.632354,0.605918,0.578976,0.552364, > 0.525820,0.499665,0.473306,0.447125,0.421027,0.394926, > 0.369017,0.343250,0.317596,0.292061,0.266805,0.241770, > 0.217121,0.192830,0.169123,0.145809,0.122984,0.100685, > 0.078967,0.057983,0.038065,0.020236,0.006364,0.000067, > 0.732957,0.720881,0.707633,0.692613,0.676200,0.657914, > 0.637626,0.615846,0.592373,0.567813,0.542877,0.518036, > 0.493386,0.468899,0.444596,0.420291,0.395908,0.371641, > 0.347423,0.323358,0.299396,0.275511,0.251830,0.228313, > 0.205146,0.182274,0.159813,0.137693,0.115983,0.094703, > 0.073974,0.053979,0.035068,0.018327,0.005580,0.000056, > 0.684126,0.674389,0.662432,0.649473,0.634783,0.618458, > 0.600173,0.580108,0.558464,0.535421,0.512051,0.488760, > 0.465717,0.442837,0.419997,0.397309,0.374393,0.351598, > 0.328943,0.306305,0.283685,0.261210,0.238875,0.216696, > 0.194737,0.173035,0.151692,0.130607,0.109864,0.089483, > 0.069628,0.050489,0.032492,0.016714,0.004940,0.000047, > 0.642684,0.634012,0.624311,0.612687,0.599751,0.584885, > 0.568165,0.549716,0.529352,0.507582,0.485673,0.463844, > 0.442031,0.420480,0.398892,0.377446,0.355834,0.334383, > 0.312903,0.291520,0.270128,0.248790,0.227632,0.206518, > 0.185624,0.164939,0.144527,0.124348,0.104415,0.084843, > 0.065754,0.047403,0.030217,0.015305,0.004396,0.000040, > 0.606784,0.599501,0.591216,0.581053,0.569262,0.555644, > 0.540298,0.522979,0.503888,0.483449,0.462725,0.441941, > 0.421308,0.400965,0.380515,0.360166,0.339667,0.319306, > 0.298896,0.278531,0.258224,0.237932,0.217694,0.197574, > 0.177567,0.157731,0.138168,0.118766,0.099587,0.080708, > 0.062310,0.044656,0.028210,0.014081,0.003935,0.000035, > 0.575250,0.569055,0.561834,0.552756,0.542358,0.530173, > 0.515880,0.499697,0.481592,0.462197,0.442485,0.422845, > 0.403137,0.383704,0.364332,0.344943,0.325467,0.306020, > 0.286496,0.267106,0.247644,0.228251,0.208881,0.189557, > 0.170377,0.151316,0.132475,0.113748,0.095226,0.076992, > 0.059214,0.042195,0.026424,0.013004,0.003541,0.000030, > 2.510860,2.424890,2.338800,2.251430,2.161030,2.070080, > 1.976500,1.884380,1.789520,1.695430,1.600900,1.507110, > 1.413980,1.321430,1.229430,1.138950,1.049380,0.961914, > 0.875863,0.792403,0.711728,0.634220,0.560088,0.490193, > 0.424574,0.363647,0.307681,0.256400,0.209654,0.167105, > 0.128273,0.092911,0.060892,0.033178,0.011804,0.000172, > 1.797890,1.740890,1.683730,1.625330,1.564220,1.502940, > 1.439440,1.375560,1.311540,1.247060,1.182260,1.117820, > 1.053940,0.990598,0.927193,0.864982,0.802868,0.741596, > 0.681232,0.621931,0.563990,0.507817,0.452982,0.400714, > 0.350809,0.303514,0.259139,0.217505,0.178849,0.142987, > 0.109689,0.079022,0.051121,0.027212,0.009208,0.000121, > 1.421860,1.380980,1.338550,1.294150,1.249120,1.201850, > 1.153940,1.105280,1.054910,1.004680,0.954677,0.905086, > 0.855141,0.806467,0.757283,0.708817,0.660443,0.613233, > 0.565935,0.519307,0.473770,0.428741,0.385002,0.342676, > 0.301797,0.262591,0.225442,0.189974,0.156683,0.125290, > 0.095921,0.068630,0.043899,0.022874,0.007411,0.000090, > 1.187690,1.155530,1.122280,1.087030,1.051020,1.013270, > 0.973852,0.933953,0.892905,0.851713,0.810455,0.769192, > 0.728379,0.687798,0.647423,0.607517,0.567475,0.528048, > 0.488850,0.450119,0.411897,0.374353,0.337466,0.301396, > 0.266522,0.232713,0.200387,0.169339,0.139725,0.111661, > 0.085284,0.060612,0.038324,0.019572,0.006096,0.000069, > 1.028240,1.002610,0.975714,0.947100,0.916605,0.885287, > 0.852214,0.817988,0.782800,0.747390,0.711939,0.676325, > 0.641102,0.606200,0.571466,0.536994,0.502546,0.468438, > 0.434547,0.400992,0.367663,0.334910,0.302578,0.270890, > 0.240102,0.210041,0.181153,0.153175,0.126406,0.100826, > 0.076605,0.054074,0.033774,0.016896,0.005055,0.000053, > 0.923239,0.901357,0.878307,0.853928,0.827994,0.800448, > 0.770833,0.740864,0.709208,0.677436,0.645505,0.613815, > 0.582211,0.551121,0.519863,0.488847,0.457919,0.427291, > 0.396728,0.366593,0.336576,0.306991,0.277689,0.248974, > 0.220914,0.193485,0.166892,0.141077,0.116297,0.092513, > 0.069960,0.049024,0.030239,0.014840,0.004277,0.000042, > 0.840109,0.821572,0.802221,0.780737,0.758007,0.733570, > 0.707391,0.680060,0.651741,0.622598,0.593525,0.564790, > 0.535923,0.507607,0.479180,0.450770,0.422703,0.394808, > 0.366922,0.339201,0.311718,0.284555,0.257727,0.231240, > 0.205291,0.179821,0.155134,0.131145,0.107877,0.085551, > 0.064423,0.044815,0.027334,0.013173,0.003668,0.000034, > 0.771898,0.756403,0.739454,0.720917,0.700462,0.678667, > 0.655360,0.630295,0.604073,0.577412,0.550894,0.524237, > 0.497613,0.471552,0.445325,0.419311,0.393354,0.367582, > 0.341816,0.316292,0.290960,0.265746,0.240859,0.216190, > 0.191974,0.168255,0.145100,0.122534,0.100715,0.079595, > 0.059666,0.041209,0.024869,0.011778,0.003176,0.000028, > 0.715432,0.701697,0.687335,0.671043,0.652774,0.633196, > 0.611618,0.588703,0.564555,0.539941,0.515046,0.490273, > 0.465666,0.441417,0.417189,0.393083,0.368842,0.344918, > 0.320853,0.296990,0.273346,0.249853,0.226554,0.203403, > 0.180727,0.158303,0.136484,0.115160,0.094448,0.074509, > 0.055551,0.038122,0.022769,0.010610,0.002775,0.000023, > 0.666960,0.655627,0.643169,0.628497,0.612317,0.594375, > 0.574670,0.553372,0.531024,0.507805,0.484742,0.461630, > 0.438595,0.415762,0.393171,0.370479,0.347803,0.325380, > 0.302957,0.280541,0.258199,0.236074,0.214141,0.192375, > 0.170880,0.149710,0.129011,0.108739,0.089084,0.069999, > 0.052020,0.035461,0.020980,0.009629,0.002450,0.000020, > 0.626089,0.616427,0.605000,0.592191,0.577470,0.561221, > 0.543098,0.523316,0.502350,0.480560,0.458749,0.436948, > 0.415413,0.393933,0.372648,0.351193,0.329853,0.308661, > 0.287323,0.266214,0.245090,0.224191,0.203371,0.182736, > 0.162324,0.142202,0.122478,0.103078,0.084237,0.066064, > 0.048902,0.033114,0.019409,0.008778,0.002175,0.000017, > 0.590664,0.582220,0.572479,0.560768,0.547507,0.532490, > 0.515657,0.497052,0.477369,0.456893,0.436199,0.415686, > 0.395192,0.374819,0.354542,0.334339,0.314120,0.293919, > 0.273897,0.253708,0.233747,0.213772,0.193977,0.174220, > 0.154737,0.135524,0.116625,0.098102,0.080014,0.062548, > 0.046083,0.031041,0.018032,0.008042,0.001942,0.000014, > 0.559453,0.552473,0.543647,0.533124,0.521070,0.507229, > 0.491452,0.474143,0.455506,0.435945,0.416256,0.396977, > 0.377326,0.358052,0.338792,0.319598,0.300262,0.281122, > 0.261883,0.242706,0.223579,0.204566,0.185577,0.166710, > 0.148088,0.129570,0.111434,0.093593,0.076187,0.059415, > 0.043600,0.029214,0.016816,0.007398,0.001744,0.000012, > 2.494430,2.407750,2.320680,2.230570,2.139410,2.045820, > 1.954300,1.859160,1.764370,1.669990,1.574490,1.481120, > 1.386830,1.294720,1.202200,1.111670,1.021960,0.934479, > 0.848525,0.765125,0.684539,0.607178,0.533502,0.463759, > 0.398700,0.338188,0.282593,0.231728,0.185473,0.143600, > 0.105860,0.072262,0.043334,0.020543,0.005838,0.000062, > 1.781830,1.724140,1.666660,1.605560,1.544230,1.481280, > 1.418440,1.353130,1.288210,1.223120,1.158330,1.093950, > 1.029490,0.965787,0.902413,0.840091,0.777784,0.716620, > 0.656180,0.597156,0.539338,0.483075,0.428748,0.376710, > 0.327154,0.280125,0.236105,0.194961,0.156820,0.121629, > 0.089495,0.060685,0.035892,0.016616,0.004508,0.000044, > 1.407480,1.365110,1.320730,1.276090,1.228930,1.181910, > 1.132930,1.082570,1.032770,0.982143,0.932084,0.882066, > 0.832465,0.783250,0.734272,0.685821,0.637505,0.589736, > 0.542652,0.496248,0.450545,0.405906,0.362451,0.320187, > 0.279650,0.240757,0.203901,0.169026,0.136216,0.105650, > 0.077490,0.052176,0.030469,0.013808,0.003599,0.000032, > 1.173640,1.140320,1.105880,1.069540,1.032320,0.993897, > 0.953913,0.913021,0.871617,0.829958,0.788517,0.747238, > 0.706258,0.666094,0.625534,0.585422,0.545512,0.506014, > 0.466972,0.428374,0.390208,0.352773,0.316057,0.280282, > 0.245648,0.212195,0.180160,0.149604,0.120645,0.093445, > 0.068267,0.045631,0.026330,0.011698,0.002940,0.000025, > 1.014900,0.987942,0.959684,0.930270,0.899044,0.866678, > 0.832642,0.797620,0.762091,0.726232,0.690987,0.655077, > 0.620091,0.585191,0.550543,0.515986,0.481432,0.447342, > 0.413584,0.380034,0.346933,0.314295,0.282265,0.250834, > 0.220232,0.190535,0.162003,0.134530,0.108400,0.083726, > 0.060878,0.040604,0.022978,0.010002,0.002422,0.000019, > 0.909560,0.887017,0.862965,0.837572,0.810211,0.781875, > 0.751746,0.720666,0.688703,0.656776,0.624847,0.593122, > 0.561398,0.530391,0.499138,0.468280,0.437229,0.406831, > 0.376391,0.346350,0.316458,0.286984,0.257979,0.229524, > 0.201720,0.174606,0.148473,0.123248,0.099111,0.076308, > 0.055184,0.036277,0.020402,0.008711,0.002038,0.000015, > 0.827323,0.808240,0.787194,0.764785,0.741102,0.715710, > 0.688724,0.660472,0.631563,0.602307,0.573344,0.544275, > 0.515643,0.487203,0.458901,0.430660,0.402505,0.374606, > 0.346855,0.319408,0.292040,0.265141,0.238553,0.212326, > 0.186666,0.161641,0.137398,0.113940,0.091458,0.070180, > 0.050481,0.032923,0.018295,0.007672,0.001739,0.000012, > 0.759389,0.742705,0.724897,0.705170,0.684041,0.661178, > 0.636646,0.610897,0.584460,0.557556,0.530715,0.504308, > 0.477779,0.451618,0.425597,0.399640,0.373702,0.348041, > 0.322356,0.297020,0.271779,0.246872,0.222152,0.197859, > 0.173964,0.150642,0.127976,0.105994,0.084903,0.064934, > 0.046463,0.030071,0.016527,0.006814,0.001499,0.000010, > 0.703019,0.688814,0.673059,0.655722,0.636535,0.615667, > 0.593245,0.569605,0.545086,0.520194,0.495411,0.470863, > 0.446201,0.421921,0.397793,0.373652,0.349528,0.325620, > 0.301845,0.278188,0.254649,0.231349,0.208318,0.185538, > 0.163153,0.141239,0.119901,0.099172,0.079281,0.060422, > 0.043027,0.027647,0.015032,0.006099,0.001305,0.000008, > 0.655549,0.643053,0.628966,0.613794,0.596440,0.577313, > 0.556691,0.534653,0.511805,0.488542,0.465448,0.442353, > 0.419418,0.396807,0.374153,0.351526,0.328993,0.306561, > 0.284259,0.262079,0.240000,0.218086,0.196361,0.174927, > 0.153790,0.133102,0.112916,0.093252,0.074414,0.056530, > 0.040069,0.025572,0.013769,0.005504,0.001148,0.000007, > 0.614886,0.603957,0.591709,0.577655,0.561834,0.544400, > 0.525193,0.504891,0.483435,0.461525,0.439729,0.418166, > 0.396478,0.375064,0.353725,0.332565,0.311308,0.290188, > 0.269125,0.248171,0.227323,0.206600,0.186070,0.165736, > 0.145686,0.125993,0.106807,0.088112,0.070133,0.053112, > 0.037464,0.023751,0.012665,0.004991,0.001015,0.000006, > 0.579569,0.570010,0.559085,0.546465,0.532125,0.516026, > 0.498089,0.478853,0.458681,0.438173,0.417404,0.396958, > 0.376499,0.356301,0.336115,0.316014,0.295930,0.275886, > 0.255871,0.236010,0.216223,0.196524,0.176991,0.157626, > 0.138515,0.119742,0.101404,0.083545,0.066339,0.050078, > 0.035167,0.022152,0.011703,0.004550,0.000904,0.000005, > 0.548753,0.540580,0.530520,0.519276,0.506004,0.491061, > 0.474436,0.456132,0.436996,0.417488,0.397755,0.378421, > 0.359085,0.339851,0.320700,0.301518,0.282357,0.263296, > 0.244273,0.225311,0.206421,0.187624,0.168959,0.150424, > 0.132145,0.114172,0.096557,0.079446,0.062967,0.047378, > 0.033126,0.020736,0.010987,0.004166,0.000809,0.000004, > 2.486670,2.401390,2.313150,2.223100,2.132210,2.040620, > 1.946240,1.852960,1.758660,1.663790,1.569300,1.475400, > 1.381630,1.288980,1.197380,1.106900,1.017420,0.930068, > 0.844430,0.760955,0.680549,0.603582,0.530075,0.460570, > 0.395641,0.335338,0.279919,0.229212,0.183101,0.141385, > 0.103826,0.070754,0.041898,0.019580,0.005430,0.000056, > 1.776150,1.719900,1.661320,1.600440,1.539110,1.477010, > 1.412790,1.348550,1.283220,1.218900,1.153830,1.088770, > 1.025080,0.961812,0.898582,0.836002,0.773940,0.712963, > 0.652650,0.593610,0.536099,0.480001,0.425916,0.373940, > 0.324454,0.277574,0.233716,0.192654,0.154624,0.119622, > 0.087676,0.059112,0.034659,0.015816,0.004189,0.000039, > 1.403340,1.361000,1.316930,1.271610,1.225380,1.177160, > 1.128240,1.078580,1.028640,0.977947,0.928063,0.878521, > 0.828550,0.779706,0.730569,0.682199,0.633948,0.586522, > 0.539401,0.493171,0.447639,0.403145,0.359673,0.317729, > 0.277222,0.238454,0.201716,0.166913,0.134251,0.103840, > 0.075864,0.050772,0.029388,0.013129,0.003341,0.000029, > 1.169950,1.136380,1.102100,1.065820,1.028190,0.989956, > 0.950189,0.908974,0.867893,0.826368,0.784940,0.743815, > 0.702872,0.662550,0.622187,0.582284,0.542431,0.503183, > 0.464112,0.425583,0.387521,0.350185,0.313576,0.277943, > 0.243398,0.210036,0.178125,0.147636,0.118821,0.091774, > 0.066794,0.044372,0.025371,0.011112,0.002727,0.000022, > 1.011370,0.984724,0.956435,0.926814,0.895443,0.862964, > 0.829274,0.794479,0.758875,0.722999,0.687335,0.652160, > 0.616882,0.582106,0.547373,0.512881,0.478596,0.444627, > 0.410813,0.377446,0.344502,0.311951,0.279953,0.248594, > 0.218124,0.188529,0.160091,0.132732,0.106696,0.082174, > 0.059499,0.039204,0.022122,0.009492,0.002245,0.000017, > 0.906689,0.883943,0.860162,0.834079,0.807067,0.778589, > 0.748446,0.717460,0.685716,0.653595,0.621726,0.590095, > 0.558631,0.527321,0.496368,0.465511,0.434680,0.404249, > 0.373821,0.343908,0.314123,0.284707,0.255819,0.227414, > 0.199673,0.172714,0.146637,0.121510,0.097519,0.074849, > 0.053902,0.035215,0.019623,0.008259,0.001888,0.000014, > 0.824084,0.805081,0.784198,0.761980,0.737877,0.712401, > 0.685185,0.657419,0.628573,0.599438,0.570456,0.541554, > 0.512854,0.484478,0.456174,0.428019,0.399930,0.372215, > 0.344523,0.317092,0.289826,0.262998,0.236396,0.210320, > 0.184715,0.159801,0.135667,0.112306,0.089938,0.068791, > 0.049274,0.031936,0.017585,0.007270,0.001610,0.000011, > 0.756740,0.739909,0.722048,0.702184,0.681164,0.658293, > 0.633526,0.608140,0.581469,0.554572,0.527904,0.501517, > 0.475145,0.448970,0.422970,0.397121,0.371248,0.345564, > 0.320059,0.294788,0.269623,0.244755,0.220148,0.195923, > 0.172134,0.148871,0.126312,0.104417,0.083450,0.063624, > 0.045337,0.029156,0.015873,0.006450,0.001387,0.000009, > 0.700680,0.686067,0.670562,0.652740,0.633579,0.612848, > 0.590476,0.566706,0.542207,0.517556,0.492696,0.468105, > 0.443562,0.419350,0.395237,0.371204,0.347185,0.323339, > 0.299543,0.276011,0.252614,0.229355,0.206336,0.183677, > 0.161386,0.139527,0.118283,0.097686,0.077903,0.059183, > 0.041956,0.026789,0.014429,0.005771,0.001207,0.000008, > 0.652991,0.640812,0.626563,0.610980,0.593505,0.574586, > 0.553939,0.531965,0.509196,0.485938,0.462796,0.439793, > 0.416876,0.394269,0.371716,0.349126,0.326720,0.304355, > 0.282067,0.259966,0.237971,0.216140,0.194490,0.173138, > 0.152089,0.131483,0.111379,0.091839,0.073081,0.055347, > 0.039055,0.024764,0.013208,0.005204,0.001061,0.000006, > 0.612493,0.601527,0.589260,0.574987,0.559293,0.541793, > 0.522714,0.502071,0.480580,0.459021,0.437202,0.415483, > 0.394052,0.372703,0.351511,0.330271,0.309080,0.288012, > 0.267062,0.246156,0.225353,0.204703,0.184235,0.163954, > 0.144013,0.124421,0.105317,0.086709,0.068861,0.051973, > 0.036498,0.022991,0.012143,0.004716,0.000938,0.000005, > 0.577224,0.567718,0.556768,0.543992,0.529614,0.513264, > 0.495593,0.476332,0.456118,0.435643,0.414958,0.394475, > 0.374159,0.353948,0.333917,0.313860,0.293761,0.273842, > 0.253865,0.234066,0.214285,0.194711,0.175202,0.155932, > 0.136903,0.118198,0.099962,0.082192,0.065107,0.048988, > 0.034247,0.021430,0.011215,0.004297,0.000835,0.000005, > 0.546580,0.538342,0.528410,0.516932,0.503593,0.488594, > 0.471844,0.453700,0.434588,0.415081,0.395496,0.376054, > 0.356702,0.337597,0.318488,0.299383,0.280225,0.261271, > 0.242274,0.223404,0.204576,0.185839,0.167243,0.148787, > 0.130586,0.112676,0.095184,0.078140,0.061781,0.046329, > 0.032244,0.020050,0.010400,0.003932,0.000747,0.000004/ DATA LQ > /0.514104,0.515832,0.517091,0.517536,0.517832,0.517221, > 0.516574,0.515224,0.513287,0.510858,0.508000,0.504917, > 0.501246,0.497403,0.492983,0.488152,0.482358,0.476122, > 0.468649,0.460339,0.450645,0.439812,0.427166,0.412836, > 0.396439,0.377753,0.356388,0.331566,0.302870,0.269780, > 0.231439,0.187668,0.138736,0.086850,0.036211,0.000604, > 0.579741,0.577476,0.575013,0.571386,0.567464,0.562828, > 0.557414,0.551410,0.544332,0.536891,0.528953,0.521049, > 0.512612,0.503967,0.494945,0.485619,0.475531,0.464958, > 0.453985,0.442176,0.429817,0.416392,0.402026,0.386267, > 0.369251,0.350413,0.329631,0.306106,0.279316,0.248539, > 0.213030,0.172492,0.127054,0.078865,0.032177,0.000513, > 0.617205,0.612889,0.607448,0.601280,0.593987,0.585725, > 0.577239,0.567645,0.557436,0.546468,0.535124,0.523690, > 0.511913,0.500192,0.488271,0.476132,0.463482,0.450653, > 0.437460,0.423920,0.409767,0.395178,0.379804,0.363620, > 0.346492,0.327958,0.307882,0.285606,0.260283,0.231441, > 0.198228,0.160224,0.117620,0.072457,0.029003,0.000445, > 0.635918,0.629589,0.622092,0.614001,0.604751,0.594717, > 0.583683,0.572041,0.559188,0.546115,0.532523,0.518899, > 0.505201,0.491782,0.477675,0.463752,0.449738,0.435623, > 0.421062,0.406484,0.391587,0.376279,0.360631,0.344349, > 0.327299,0.309280,0.289970,0.268596,0.244665,0.217365, > 0.185991,0.150104,0.109826,0.067166,0.026423,0.000392, > 0.642211,0.634584,0.625818,0.616705,0.605951,0.594718, > 0.582351,0.569037,0.554986,0.540276,0.525303,0.510521, > 0.495596,0.480637,0.465827,0.451058,0.435872,0.421066, > 0.405880,0.390796,0.375397,0.359903,0.344185,0.327987, > 0.311213,0.293614,0.274949,0.254454,0.231575,0.205516, > 0.175652,0.141477,0.103138,0.062618,0.024219,0.000348, > 0.639180,0.631092,0.621893,0.612137,0.600935,0.589006, > 0.576122,0.561754,0.547048,0.531226,0.515512,0.499878, > 0.484089,0.468662,0.453258,0.437735,0.422292,0.407026, > 0.391544,0.376258,0.360735,0.345208,0.329470,0.313579, > 0.297158,0.280066,0.261996,0.242200,0.220305,0.195335, > 0.166765,0.134023,0.097367,0.058717,0.022356,0.000312, > 0.631156,0.622788,0.613580,0.603527,0.592277,0.579779, > 0.566484,0.551813,0.536423,0.520508,0.504251,0.488021, > 0.472093,0.456329,0.440551,0.424908,0.409169,0.393745, > 0.378229,0.362803,0.347473,0.331992,0.316442,0.300771, > 0.284762,0.268174,0.250682,0.231645,0.210500,0.186498, > 0.158976,0.127569,0.092372,0.055348,0.020770,0.000282, > 0.619698,0.611519,0.602336,0.591943,0.580814,0.568514, > 0.554982,0.540509,0.524970,0.508623,0.492412,0.476046, > 0.459935,0.444036,0.428166,0.412462,0.396815,0.381305, > 0.365892,0.350546,0.335395,0.320121,0.304877,0.289437, > 0.273841,0.257707,0.240774,0.222305,0.201894,0.178726, > 0.152209,0.121898,0.087984,0.052399,0.019398,0.000257, > 0.606343,0.598499,0.589365,0.579821,0.568434,0.556327, > 0.542872,0.528454,0.512900,0.496651,0.480390,0.464149, > 0.447972,0.432310,0.416261,0.400748,0.385116,0.369727, > 0.354542,0.339434,0.324280,0.309328,0.294433,0.279362, > 0.264086,0.248363,0.231853,0.214025,0.194261,0.171833, > 0.146197,0.116876,0.084113,0.049805,0.018206,0.000236, > 0.591896,0.584256,0.575797,0.566248,0.555494,0.543797, > 0.530664,0.516396,0.501105,0.484951,0.468541,0.452368, > 0.436358,0.420764,0.405075,0.389617,0.374238,0.359020, > 0.343972,0.329043,0.314268,0.299607,0.284884,0.270198, > 0.255337,0.239949,0.223985,0.206644,0.187466,0.165701, > 0.140843,0.112412,0.080664,0.047512,0.017165,0.000218, > 0.576697,0.569585,0.561998,0.552972,0.542790,0.531174, > 0.518335,0.504287,0.489160,0.473285,0.457139,0.441230, > 0.425520,0.409903,0.394534,0.379177,0.364141,0.349046, > 0.334214,0.319602,0.305086,0.290690,0.276241,0.261890, > 0.247381,0.232411,0.216817,0.199969,0.181340,0.160182, > 0.135997,0.108375,0.077551,0.045446,0.016236,0.000202, > 0.561779,0.555491,0.548244,0.539719,0.529917,0.518731, > 0.506334,0.492791,0.477798,0.462049,0.446208,0.430534, > 0.414938,0.399619,0.384500,0.369572,0.354611,0.339912, > 0.325296,0.310898,0.296649,0.282465,0.268439,0.254375, > 0.240158,0.225578,0.210376,0.193910,0.175777,0.155185, > 0.131618,0.104727,0.074745,0.043586,0.015406,0.000188, > 0.547406,0.541522,0.534772,0.526637,0.517320,0.506665, > 0.494929,0.481541,0.466958,0.451581,0.435977,0.420382, > 0.405163,0.390039,0.375242,0.360562,0.345699,0.331352, > 0.316934,0.302875,0.288897,0.275075,0.261211,0.247494, > 0.233546,0.219335,0.204477,0.188465,0.170694,0.150613, > 0.127623,0.101395,0.072188,0.041898,0.014659,0.000176, > 0.510949,0.512461,0.513125,0.513775,0.513627,0.512953, > 0.511676,0.510153,0.508026,0.505586,0.502488,0.499514, > 0.495918,0.491777,0.487221,0.482318,0.476322,0.469865, > 0.462254,0.453539,0.443579,0.431957,0.418405,0.402714, > 0.384513,0.363211,0.338441,0.309525,0.276023,0.237418, > 0.194095,0.146991,0.098692,0.053921,0.018675,0.000251, > 0.573840,0.571539,0.568660,0.564636,0.560125,0.555429, > 0.549250,0.542654,0.535827,0.528339,0.520207,0.512007, > 0.503505,0.494854,0.485721,0.476289,0.466253,0.455649, > 0.444534,0.432471,0.419685,0.405709,0.390526,0.373656, > 0.355038,0.334208,0.310578,0.283443,0.252376,0.217005, > 0.177144,0.133976,0.089700,0.048645,0.016531,0.000213, > 0.610274,0.604837,0.598760,0.591967,0.584049,0.575825, > 0.566634,0.556589,0.545957,0.534847,0.523533,0.511998, > 0.500328,0.488702,0.476646,0.464532,0.451851,0.438900, > 0.425698,0.412022,0.397498,0.382650,0.366763,0.349574, > 0.330976,0.310680,0.288132,0.262614,0.233680,0.200673, > 0.163728,0.123628,0.082475,0.044455,0.014854,0.000185, > 0.627192,0.620322,0.611968,0.603362,0.593426,0.582683, > 0.571079,0.558798,0.545824,0.532347,0.519069,0.505318, > 0.491597,0.478201,0.464023,0.450314,0.436183,0.422097, > 0.407576,0.392963,0.377856,0.362363,0.346096,0.329071, > 0.310942,0.291304,0.269699,0.245618,0.218210,0.187279, > 0.152695,0.115069,0.076563,0.041003,0.013496,0.000163, > 0.631625,0.623836,0.614536,0.604467,0.593362,0.581189, > 0.568126,0.554473,0.540005,0.525342,0.510320,0.495556, > 0.480543,0.465856,0.450901,0.436047,0.421171,0.406166, > 0.391154,0.375993,0.360616,0.344868,0.328695,0.311826, > 0.294125,0.275057,0.254485,0.231491,0.205428,0.176164, > 0.143362,0.107930,0.071534,0.038049,0.012338,0.000145, > 0.628674,0.619289,0.609887,0.598892,0.587603,0.574461, > 0.560634,0.545957,0.530648,0.515019,0.499229,0.483499, > 0.467857,0.452587,0.437187,0.421817,0.406489,0.391076, > 0.375987,0.360550,0.345003,0.329399,0.313389,0.296812, > 0.279488,0.261191,0.241276,0.219209,0.194576,0.166551, > 0.135508,0.101696,0.067217,0.035540,0.011363,0.000130, > 0.619655,0.610842,0.600741,0.589506,0.577444,0.564098, > 0.550163,0.534861,0.519180,0.503427,0.486955,0.471011, > 0.455040,0.439305,0.423732,0.408132,0.392629,0.377169, > 0.361831,0.346498,0.331073,0.315724,0.299859,0.283638, > 0.266810,0.249063,0.229985,0.208716,0.185019,0.158309, > 0.128580,0.096355,0.063475,0.033376,0.010535,0.000117, > 0.607474,0.598863,0.589070,0.578071,0.565959,0.552353, > 0.538238,0.523141,0.507028,0.491211,0.474478,0.458296, > 0.442260,0.426619,0.410731,0.395250,0.379629,0.364356, > 0.348971,0.333841,0.318608,0.303265,0.287813,0.272025, > 0.255702,0.238406,0.219910,0.199549,0.176755,0.151092, > 0.122529,0.091666,0.060228,0.031496,0.009822,0.000107, > 0.593760,0.585285,0.576038,0.564716,0.553007,0.539720, > 0.525715,0.510544,0.494716,0.478059,0.461888,0.445882, > 0.429808,0.413922,0.398545,0.382937,0.367557,0.352300, > 0.337238,0.322197,0.307168,0.292182,0.277095,0.261674, > 0.245761,0.228939,0.211113,0.191442,0.169413,0.144749, > 0.117225,0.087581,0.057361,0.029844,0.009202,0.000098, > 0.579324,0.570883,0.561883,0.551297,0.539642,0.526840, > 0.512914,0.497903,0.481998,0.465838,0.449788,0.433693, > 0.417703,0.402286,0.386792,0.371437,0.356215,0.341283, > 0.326313,0.311637,0.296947,0.282140,0.267321,0.252270, > 0.236761,0.220580,0.203236,0.184230,0.162954,0.139137, > 0.112557,0.083990,0.054824,0.028391,0.008663,0.000090, > 0.564315,0.556562,0.547782,0.537905,0.526572,0.514011, > 0.500033,0.485378,0.470013,0.453915,0.438022,0.422153, > 0.406530,0.391112,0.375843,0.360792,0.345939,0.331089, > 0.316383,0.301913,0.287442,0.273075,0.258564,0.243883, > 0.228827,0.213013,0.196181,0.177712,0.157167,0.134036, > 0.108376,0.080670,0.052555,0.027085,0.008183,0.000084, > 0.549474,0.542540,0.534001,0.524222,0.513538,0.501429, > 0.487998,0.473665,0.458530,0.442760,0.426794,0.411187, > 0.395839,0.380807,0.365652,0.350917,0.336154,0.321565, > 0.307328,0.293072,0.278886,0.264801,0.250668,0.236339, > 0.221517,0.206147,0.189835,0.171902,0.151860,0.129450, > 0.104545,0.077707,0.050505,0.025913,0.007753,0.000078, > 0.535066,0.528539,0.520546,0.511350,0.500893,0.489256, > 0.476407,0.462092,0.447135,0.431648,0.416302,0.400932, > 0.385885,0.370871,0.356128,0.341562,0.327354,0.312933, > 0.298890,0.284864,0.270948,0.257177,0.243326,0.229321, > 0.214967,0.199986,0.183994,0.166574,0.147078,0.125313, > 0.101089,0.075065,0.048640,0.024855,0.007368,0.000073, > 0.509409,0.510722,0.511674,0.511738,0.511286,0.510766, > 0.509198,0.507564,0.505152,0.502704,0.499650,0.496452, > 0.492756,0.488747,0.483986,0.478817,0.472835,0.466164, > 0.458412,0.449362,0.438892,0.426796,0.412430,0.395391, > 0.375497,0.352110,0.324651,0.292275,0.254919,0.212461, > 0.165886,0.117521,0.071431,0.033526,0.009193,0.000091, > 0.570888,0.568118,0.564921,0.561074,0.555970,0.550478, > 0.544438,0.537826,0.530394,0.523005,0.514685,0.506342, > 0.497778,0.489102,0.479852,0.470348,0.460202,0.449558, > 0.438121,0.426046,0.412891,0.398444,0.382493,0.364765, > 0.344689,0.321925,0.295886,0.265915,0.231610,0.192817, > 0.150406,0.106352,0.064426,0.030028,0.008099,0.000077, > 0.605452,0.599834,0.593296,0.586627,0.578039,0.569200, > 0.559754,0.549624,0.538672,0.527535,0.516054,0.504396, > 0.492657,0.481104,0.468974,0.456682,0.443949,0.431100, > 0.417683,0.403750,0.389161,0.373852,0.357151,0.339250, > 0.319633,0.297684,0.273111,0.245042,0.213167,0.177272, > 0.138115,0.097498,0.058870,0.027271,0.007249,0.000067, > 0.621534,0.614157,0.605835,0.596272,0.585946,0.574335, > 0.562822,0.550233,0.536841,0.523598,0.509920,0.496293, > 0.482476,0.468781,0.455150,0.441275,0.427053,0.412946, > 0.398426,0.383702,0.368436,0.352438,0.335827,0.317985, > 0.298878,0.277758,0.254411,0.228040,0.198146,0.164659, > 0.128112,0.090261,0.054331,0.025017,0.006563,0.000059, > 0.625646,0.616575,0.607120,0.596395,0.584780,0.572447, > 0.558731,0.544391,0.529948,0.515041,0.499832,0.485344, > 0.470145,0.455555,0.440783,0.425851,0.410922,0.396061, > 0.381038,0.365825,0.350234,0.334195,0.317648,0.300202, > 0.281570,0.261329,0.238985,0.213898,0.185675,0.154079, > 0.119726,0.084159,0.050488,0.023108,0.005981,0.000052, > 0.621816,0.611958,0.601503,0.590515,0.577713,0.564667, > 0.550311,0.535181,0.519529,0.503933,0.488235,0.472641, > 0.456837,0.441572,0.426210,0.410935,0.395680,0.380387, > 0.365098,0.349710,0.334114,0.318243,0.301798,0.284744, > 0.266642,0.247149,0.225778,0.201848,0.175065,0.145114, > 0.112575,0.078976,0.047163,0.021485,0.005492,0.000047, > 0.612456,0.602615,0.592112,0.580629,0.567547,0.554000, > 0.539266,0.523505,0.507494,0.491317,0.475170,0.459335, > 0.443356,0.427588,0.412089,0.396613,0.381018,0.365836, > 0.350418,0.335168,0.319712,0.303995,0.287942,0.271269, > 0.253737,0.234846,0.214374,0.191514,0.165911,0.137420, > 0.106436,0.074522,0.044401,0.020096,0.005079,0.000042, > 0.599906,0.590658,0.579871,0.568198,0.555552,0.541752, > 0.526573,0.511088,0.494690,0.478321,0.462067,0.446064, > 0.430078,0.414174,0.398700,0.383180,0.367731,0.352392, > 0.337125,0.322012,0.306729,0.291406,0.275582,0.259369, > 0.242384,0.224158,0.204420,0.182500,0.157984,0.130669, > 0.101087,0.070603,0.041960,0.018893,0.004724,0.000039, > 0.586385,0.576832,0.566175,0.554843,0.542371,0.528570, > 0.513481,0.497995,0.481584,0.465454,0.449166,0.432955, > 0.417165,0.401569,0.386012,0.370517,0.355274,0.340119, > 0.324992,0.310107,0.295066,0.280058,0.264679,0.248876, > 0.232336,0.214712,0.195681,0.174542,0.150960,0.124776, > 0.096383,0.067222,0.039821,0.017842,0.004416,0.000035, > 0.571562,0.562230,0.553483,0.541106,0.528820,0.515154, > 0.500724,0.484959,0.468895,0.452796,0.436478,0.420645, > 0.404877,0.389321,0.373971,0.358762,0.343687,0.328774, > 0.313968,0.299216,0.284619,0.269828,0.254780,0.239398, > 0.223345,0.206321,0.187886,0.167537,0.144752,0.119523, > 0.092229,0.064212,0.037929,0.016920,0.004149,0.000033, > 0.556430,0.547716,0.538118,0.527507,0.515330,0.502167, > 0.487523,0.472378,0.456637,0.440513,0.424536,0.408900, > 0.393308,0.378059,0.362874,0.347899,0.333001,0.318366, > 0.303866,0.289409,0.275001,0.260570,0.245889,0.230899, > 0.215327,0.198784,0.180938,0.161191,0.139240,0.114867, > 0.088515,0.061527,0.036242,0.016094,0.003912,0.000030, > 0.541088,0.533114,0.524472,0.513843,0.502197,0.489243, > 0.475392,0.460194,0.444773,0.428968,0.413337,0.397800, > 0.382488,0.367387,0.352568,0.337799,0.323189,0.308771, > 0.294545,0.280383,0.266247,0.252107,0.237863,0.223231, > 0.208050,0.191998,0.174706,0.155516,0.134237,0.110630, > 0.085174,0.059101,0.034727,0.015356,0.003700,0.000028, > 0.526746,0.519306,0.510794,0.500680,0.489449,0.477026, > 0.463177,0.448622,0.433350,0.417867,0.402498,0.387171, > 0.372235,0.357466,0.342866,0.328463,0.314082,0.299986, > 0.285991,0.272106,0.258305,0.244422,0.230511,0.216246, > 0.201466,0.185822,0.169011,0.150402,0.129717,0.106819, > 0.082146,0.056897,0.033348,0.014690,0.003511,0.000026, > 0.508294,0.509537,0.510232,0.510364,0.510028,0.509066, > 0.507574,0.506096,0.503764,0.501269,0.498261,0.494940, > 0.491357,0.487139,0.482585,0.477518,0.471487,0.464685, > 0.456954,0.448138,0.437547,0.425359,0.410953,0.393893, > 0.373877,0.350282,0.322653,0.290145,0.252556,0.209889, > 0.163221,0.114885,0.069163,0.031966,0.008546,0.000081, > 0.569358,0.566803,0.562983,0.559001,0.554723,0.548955, > 0.542742,0.536150,0.528626,0.520981,0.513074,0.504898, > 0.496364,0.487653,0.478428,0.468815,0.458593,0.447996, > 0.436563,0.424615,0.411428,0.396895,0.381014,0.363131, > 0.342953,0.320173,0.294051,0.263852,0.229272,0.190367, > 0.147875,0.103902,0.062322,0.028615,0.007526,0.000069, > 0.603670,0.598438,0.591834,0.584727,0.576222,0.567680, > 0.558196,0.547591,0.536787,0.525654,0.514148,0.502844, > 0.491102,0.479283,0.467390,0.455046,0.442422,0.429358, > 0.416062,0.402212,0.387705,0.372222,0.355668,0.337678, > 0.318005,0.295922,0.271145,0.242997,0.210936,0.174968, > 0.135724,0.095201,0.057278,0.025971,0.006733,0.000060, > 0.619938,0.612660,0.603237,0.594399,0.584220,0.572945, > 0.560906,0.547974,0.535000,0.521640,0.508074,0.494503, > 0.480788,0.467052,0.453309,0.439430,0.425326,0.411154, > 0.396834,0.381957,0.366731,0.350903,0.334155,0.316435, > 0.297197,0.275976,0.252539,0.225982,0.196020,0.162415, > 0.125838,0.088092,0.052512,0.023816,0.006094,0.000053, > 0.623752,0.615162,0.604853,0.594199,0.582617,0.570319, > 0.556672,0.542512,0.527993,0.513108,0.498185,0.483234, > 0.468378,0.453629,0.438899,0.424139,0.409158,0.394331, > 0.379329,0.364110,0.348592,0.332662,0.316045,0.298521, > 0.279860,0.259548,0.237136,0.211942,0.183600,0.151947, > 0.117534,0.082104,0.048767,0.021982,0.005551,0.000047, > 0.619727,0.610281,0.600069,0.588100,0.575830,0.562394, > 0.548038,0.533038,0.517667,0.501730,0.486231,0.470564, > 0.454906,0.439606,0.424254,0.409150,0.393803,0.378587, > 0.363406,0.347957,0.332413,0.316594,0.300201,0.283135, > 0.264987,0.245411,0.223977,0.199987,0.173072,0.143020, > 0.110486,0.077018,0.045586,0.020431,0.005096,0.000042, > 0.610495,0.600588,0.589924,0.578181,0.565270,0.551795, > 0.537196,0.521624,0.505437,0.489150,0.472983,0.457321, > 0.441422,0.425796,0.410213,0.394710,0.379352,0.363971, > 0.348802,0.333422,0.317927,0.302300,0.286286,0.269619, > 0.252021,0.233168,0.212605,0.189631,0.163957,0.135368, > 0.104416,0.072632,0.042854,0.019101,0.004711,0.000038, > 0.597856,0.588522,0.577800,0.566289,0.553520,0.539436, > 0.524435,0.508812,0.492453,0.476231,0.459950,0.443865, > 0.427942,0.412263,0.396763,0.381179,0.365851,0.350734, > 0.335459,0.320164,0.305048,0.289712,0.274016,0.257736, > 0.240735,0.222511,0.202688,0.180671,0.156046,0.128701, > 0.099140,0.068812,0.040482,0.017948,0.004381,0.000035, > 0.584084,0.574752,0.564446,0.552705,0.540225,0.526223, > 0.511434,0.495863,0.479694,0.463424,0.447226,0.431058, > 0.415167,0.399560,0.384026,0.368648,0.353514,0.338375, > 0.323406,0.308441,0.293400,0.278400,0.263020,0.247257, > 0.230713,0.213051,0.193943,0.172756,0.149112,0.122832, > 0.094502,0.065475,0.038404,0.016946,0.004094,0.000032, > 0.569517,0.560385,0.550236,0.539113,0.526421,0.513222, > 0.498378,0.482790,0.466696,0.450744,0.434441,0.418603, > 0.402901,0.387328,0.372129,0.356950,0.341855,0.326972, > 0.312231,0.297539,0.282901,0.268148,0.253147,0.237801, > 0.221735,0.204645,0.186195,0.165769,0.142932,0.117672, > 0.090399,0.062530,0.036573,0.016065,0.003846,0.000029, > 0.554361,0.545724,0.536138,0.525127,0.513348,0.500011, > 0.485538,0.470383,0.454424,0.438520,0.422677,0.406888, > 0.391373,0.376154,0.360890,0.345988,0.331296,0.316666, > 0.302083,0.287727,0.273311,0.258904,0.244293,0.229265, > 0.213676,0.197142,0.179286,0.159481,0.137449,0.113024, > 0.086731,0.059884,0.034929,0.015279,0.003625,0.000027, > 0.539566,0.531375,0.522021,0.511690,0.499982,0.487119, > 0.474212,0.458150,0.442576,0.426822,0.411177,0.395600, > 0.380520,0.365481,0.350710,0.336033,0.321423,0.307074, > 0.292885,0.278746,0.264660,0.250555,0.236263,0.221649, > 0.206449,0.190350,0.173011,0.154255,0.132811,0.108865, > 0.083444,0.057506,0.033461,0.014573,0.003428,0.000025, > 0.524865,0.517319,0.508721,0.498692,0.487348,0.474698, > 0.461245,0.446504,0.431299,0.416072,0.400507,0.385278, > 0.370252,0.355605,0.341042,0.326661,0.312354,0.298301, > 0.284289,0.270473,0.256711,0.242903,0.228985,0.214684, > 0.199885,0.184233,0.167367,0.148740,0.128003,0.105093, > 0.080442,0.055358,0.032128,0.013936,0.003252,0.000024/ C IF (X.GE.XMAXX.OR.Q2.LT.Q2MIN.OR.Q2.GT.Q2MAX) THEN C WRITE(*,*) 'X or Q2 out of range!!' C WRITE(*,*) 'Valid range is:' C WRITE(*,*) '1E-05 < X < 0.95 and' C WRITE(*,*) '4 < Q2 < 520 GeV**2' C RETURN C END IF IF (X.GE.XMAX.OR.B.GT.BMAX) THEN RES(1) = 1.0 RES(2) = 1.0 RETURN END IF IF(Q2.LT.Q2MIN) THEN Q2 = Q2MIN ELSE IF(Q2.GT.Q2MAX) THEN Q2 = Q2MAX ENDIF CALL GGINTER(INUCL,X,Q2,SHAD) DO I=1,31 TATMP(I) = TA(I,INUCL) ENDDO CALL SPLINE(TMAX,IMPAR,TATMP,C,D,E) TAF = SEVAL(TMAX,B,IMPAR,TATMP,C,D,E) DO IK=1,2 RES(IK) = 1.0/(1.0 + (ANUCL(INUCL)-1.0)*SHAD(IK)*TAF) ENDDO RETURN END SUBROUTINE GGINTER(INUCL,X,Q2,SHAD) IMPLICIT DOUBLE PRECISION(A-H,L-Z) IMPLICIT INTEGER(I-K) DIMENSION SHAD(2) DIMENSION G(36,13,4),LQ(36,13,4) DIMENSION XB(36),Q2V(13) DIMENSION GQTMP(13),LQQTMP(13) DIMENSION GXTMP(36),LQXTMP(36) DIMENSION C(100),D(100),E(100) PARAMETER(INMAX=13,IMMAX=36) COMMON/GG07/ XB,Q2V,G,LQ DO K=1,36 DO J=1,13 GQTMP(J) = G(K,J,INUCL) LQQTMP(J) = LQ(K,J,INUCL) ENDDO CALL SPLINE(INMAX,Q2V,GQTMP,C,D,E) GXTMP(K) = SEVAL(INMAX,Q2,Q2V,GQTMP,C,D,E) CALL SPLINE(INMAX,Q2V,LQQTMP,C,D,E) LQXTMP(K) = SEVAL(INMAX,Q2,Q2V,LQQTMP,C,D,E) ENDDO CALL SPLINE(IMMAX,XB,GXTMP,C,D,E) SHAD(1) = SEVAL(IMMAX,X,XB,GXTMP,C,D,E) CALL SPLINE(IMMAX,XB,LQXTMP,C,D,E) SHAD(2) = SEVAL(IMMAX,X,XB,LQXTMP,C,D,E) RETURN END C --------------------------------------------------------------------- SUBROUTINE SPLINE(N,X,Y,B,C,D) C --------------------------------------------------------------------- c*************************************************************************** C CALCULATE THE COEFFICIENTS B,C,D IN A CUBIC SPLINE INTERPOLATION. C INTERPOLATION SUBROUTINES ARE TAKEN FROM C G.E. FORSYTHE, M.A. MALCOLM AND C.B. MOLER, C COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS (PRENTICE-HALL, 1977). IMPLICIT double precision (A-H,O-Z) DIMENSION X(100),Y(100),B(100),C(100),D(100) NM1=N-1 IF(N.LT.2) RETURN IF(N.LT.3) GO TO 250 D(1)=X(2)-X(1) C(2)=(Y(2)-Y(1))/D(1) DO 210 I=2,NM1 D(I)=X(I+1)-X(I) B(I)=2.0D0*(D(I-1)+D(I)) C(I+1)=(Y(I+1)-Y(I))/D(I) C(I)=C(I+1)-C(I) 210 CONTINUE B(1)=-D(1) B(N)=-D(N-1) C(1)=0.0D0 C(N)=0.0D0 IF(N.EQ.3) GO TO 215 C(1)=C(3)/(X(4)-X(2))-C(2)/(X(3)-X(1)) C(N)=C(N-1)/(X(N)-X(N-2))-C(N-2)/(X(N-1)-X(N-3)) C(1)=C(1)*D(1)**2.0D0/(X(4)-X(1)) C(N)=-C(N)*D(N-1)**2.0D0/(X(N)-X(N-3)) 215 CONTINUE DO 220 I=2,N T=D(I-1)/B(I-1) B(I)=B(I)-T*D(I-1) C(I)=C(I)-T*C(I-1) 220 CONTINUE C(N)=C(N)/B(N) DO 230 IB=1,NM1 I=N-IB C(I)=(C(I)-D(I)*C(I+1))/B(I) 230 CONTINUE B(N)=(Y(N)-Y(NM1))/D(NM1)+D(NM1)*(C(NM1)+2.0D0*C(N)) DO 240 I=1,NM1 B(I)=(Y(I+1)-Y(I))/D(I)-D(I)*(C(I+1)+2.0D0*C(I)) D(I)=(C(I+1)-C(I))/D(I) C(I)=3.0D0*C(I) 240 CONTINUE C(N)=3.0D0*C(N) D(N)=D(N-1) RETURN 250 CONTINUE B(1)=(Y(2)-Y(1))/(X(2)-X(1)) C(1)=0.0D0 D(1)=0.0D0 B(2)=B(1) C(2)=0.0D0 D(2)=0.0D0 RETURN END c c*************************************************************************** C --------------------------------------------------------------------- double precision FUNCTION SEVAL(N,XX,X,Y,B,C,D) C --------------------------------------------------------------------- c*************************************************************************** C CALCULATE THE DISTRIBUTION AT XX BY CUBIC SPLINE INTERPOLATION. implicit double precision(A-H,O-Z) DIMENSION X(100),Y(100),B(100),C(100),D(100) DATA I/1/ IF(I.GE.N) I=1 IF(XX.LT.X(I)) GO TO 310 IF(XX.LE.X(I+1)) GO TO 330 310 CONTINUE I=1 J=N+1 320 CONTINUE K=(I+J)/2 IF(XX.LT.X(K)) J=K IF(XX.GE.X(K)) I=K IF(J.GT.I+1) GO TO 320 330 CONTINUE DX=XX-X(I) SEVAL=Y(I)+DX*(B(I)+DX*(C(I)+DX*D(I))) RETURN END c c***************************************************************************