subroutine wtocalljs implicit real*8(a-h,o-z) * parameter (nout=99) parameter(nmxhep=2000) parameter (itmxmx=100000) * common/wtgnum/itmx common/wttopt/ios,iosf common/hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep), # jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep) * dimension aidhep(6) dimension aphep(5,6) dimension amps(itmxmx),amps1(itmxmx),amps2(itmxmx) * if(iosf.gt.0) then do i=1,6 isthep(i)= 1 enddo else do i=1,4 isthep(i)= 1 enddo endif open(nout,file='mom',status='old') * do i=1,itmx do i=1,4 nevhep= i read(nout,fmt=*) nhep read(nout,fmt=*) amps(i),amps1(i),amps2(i) if(iosf.gt.0) then read(nout,fmt=*) (aidhep(j),j=1,6) else read(nout,fmt=*) (aidhep(j),j=1,4) endif read(nout,fmt=*) (aphep(j,1),j=1,5) read(nout,fmt=*) (aphep(j,2),j=1,5) read(nout,fmt=*) (aphep(j,3),j=1,5) read(nout,fmt=*) (aphep(j,4),j=1,5) if(iosf.gt.0) then read(nout,fmt=*) (aphep(j,5),j=1,5) read(nout,fmt=*) (aphep(j,6),j=1,5) endif do j=1,4 idhep(j)= aidhep(j) enddo if(iosf.gt.0) then idhep(5)= aidhep(5) idhep(6)= aidhep(6) endif do j=1,5 do l=1,4 phep(j,l)= aphep(j,l) enddo if(iosf.gt.0) then phep(j,5)= aphep(j,5) phep(j,6)= aphep(j,6) endif enddo * istrat= 0 atotsq= amps(i) a1sq= amps1(i) a2sq= amps2(i) call lu4frm(atotsq,a1sq,a2sq,istrat,irad,itau,ierr) call lulist(1) enddo close(nout) * return end * *--------------------------------------------------------------------------- * *-- Author : * The following subroutine illustrates how to write a generic * interface between electroweak generators and QCD generators * for LEP 2 applications. * * Author of this routine: T.Sjostrand * * In this particular case, an electroweak generator is supposed to * have produced two fermions, two antifermions and an arbitrary * number of photons. These particles are stored in the HEPEVT * common block. The allowed order is specified by a standard. * In brief, the final fermions should appear in the order * fermion (1) - antifermion (2) - fermion (3) - antifermion (4). * The flavour pairs should be arranged so that, if possible, the * first two could come from a W+ and the second two from a W-; * else each pair should have flavours consistent with a Z0. * The subroutine LU4FRM is supposed to read the configuration, * and call JETSET to do parton showers and fragmentation. * * Since the colour flow need not be unique, three real numbers * should be given when LU4FRM is called: * ATOTSQ = total squared amplitude for the event, irrespective of * colour flow; * A1SQ = squared amplitude for the configuration with 1 + 2 and * 3 + 4 as the two colour singlets; and * A2SQ = squared amplitude for the configuration with 1 + 4 and * 3 + 2 as the two colour singlets. * * The choice of strategy is determined by an integer input: * ISTRAT = 0 : pick configurations according to A1SQ : A2SQ; * = 1 : assign interference to maximize 1 + 2 and 3 + 4; or * = 2 : assign interference to maximize 1 + 4 and 3 + 2. * * Final-state QED radiation may be allowed or inhibited: * IRAD = 0 : no final-state photon radiation. * = 1 : photon radiation inside each final fermion pair. * * tau lepton decay may be handled by QCD generator or not. * ITAU = 0 : taus are considered stable by QCD generator. * = 1 : taus are allowed to decay by QCD generator. * * IERR is an error flag, used both as input and output. * At input, 0 means leave routine in case of error, * nonzero means stop program execution. * At output, 0 means acceptable fermion configuration, * nonzero means treatment aborted for some reason. * It is up to the writer of the main program to pick error strategy, * i.e. to let the program crash or try to fix errors. subroutine lu4frm(atotsq,a1sq,a2sq,istrat,irad,itau,ierr) * real*8 atotsq,a1sq,a2sq * common/lujets/n,k(4000,5),p(4000,5),v(4000,5) common/ludat1/mstu(200),paru(200),mstj(200),parj(200) * dimension ijoin(2) * external rlu * * call luhepc to convert from hepevt to lujets common. * inerr= 0 mstu(28)= 0 call luhepc(2) if(mstu(28).eq.8) inerr= 1 * * loop through entries and pick up all final fermions/antifermions. * i1= 0 i2= 0 i3= 0 i4= 0 do 100 i=1,n if(k(i,1).le.0.or.k(i,1).gt.10) goto 100 kfa= iabs(k(i,2)) if((kfa.ge.1.and.kfa.le.6).or.(kfa.ge.11.and.kfa.le.16)) then if(k(i,2).gt.0) then if(i1.eq.0) then i1= i elseif(i3.eq.0) then i3= i else inerr= 2 call luerrm(6,'(lu4frm:) more than two fermions') endif else if(i2.eq.0) then i2= i elseif(i4.eq.0) then i4= i else inerr= 3 call luerrm(6,'(lu4frm:) more than two antifermions') endif endif endif 100 continue * * check that event is arranged according to conventions. * if(i1.eq.0.or.i2.eq.0.or.i3.eq.0.or.i4.eq.0) then inerr= 4 call luerrm(6,'(lu4frm:) event contains too few fermions') endif if(i2.lt.i1.or.i3.lt.i2.or.i4.lt.i3) then inerr= 5 call luerrm(6,'(lu4frm:) fermions arranged in wrong order') endif * * check which fermion pairs are quarks and which leptons. * if(iabs(k(i1,2)).lt.10.and.iabs(k(i2,2)).lt.10) then iql12= 1 elseif(iabs(k(i1,2)).gt.10.and.iabs(k(i2,2)).gt.10) then iql12= 2 else inerr= 6 call luerrm(6,'(lu4frm:) first fermion pair inconsistent') endif if(iabs(k(i3,2)).lt.10.and.iabs(k(i4,2)).lt.10) then iql34= 1 elseif(iabs(k(i3,2)).gt.10.and.iabs(k(i4,2)).gt.10) then iql34= 2 else inerr= 7 call luerrm(6,'(lu4frm:) second fermion pair inconsistent') endif * * return or stop program in case of problems. * if(inerr.eq.0) then ierr= 0 elseif(ierr.eq.0) then ierr= inerr return else write(6,*) ' error: listing of faulty event follows:' call lulist(2) write(6,*) ' fermions found in lines ',i1,i2,i3,i4 write(6,*) ' error type in event above is ',inerr write(6,*) ' program execution will be stopped now' write(6,*) ' since main program does not correct errors!' stop endif * * decide whether to allow or not photon radiation in showers. * mstj(41)= 2 if(irad.eq.0) mstj(41)= 1 * * decide on colour pairing. * if(iql12.eq.2.and.iql34.eq.2) then npair= 0 elseif(iql12.eq.1.and.iql34.eq.2) then npair= 1 ip1= i1 ip2= i2 elseif(iql12.eq.2.and.iql34.eq.1) then npair= 1 ip1= i3 ip2= i4 else npair= 2 ip1= i1 ip3= i3 r1sq= real(a1sq) r2sq= real(a2sq) delta= real(atotsq-a1sq-a2sq) if(istrat.eq.1) then if(delta.gt.0.) r1sq= r1sq+delta if(delta.lt.0.) r2sq= max(0.,r2sq+delta) elseif(istrat.eq.2) then if(delta.gt.0.) r2sq= r2sq+delta if(delta.lt.0.) r1sq= max(0.,r1sq+delta) endif if(r2sq.lt.rlu(0)*(r1sq+r2sq)) then ip2= i2 ip4= i4 else ip2= i4 ip4= i2 endif endif * * do colour joining and parton showers. * if(npair.ge.1) then ijoin(1)= ip1 ijoin(2)= ip2 call lujoin(2,ijoin) pm12s= (p(ip1,4)+p(ip2,4))**2-(p(ip1,1)+p(ip2,1))**2- # (p(ip1,2)+p(ip2,2))**2-(p(ip1,3)+p(ip2,3))**2 call lushow(ip1,ip2,sqrt(max(0.,pm12s))) endif if(npair.eq.2) then ijoin(1)= ip3 ijoin(2)= ip4 call lujoin(2,ijoin) pm34s= (p(ip3,4)+p(ip4,4))**2-(p(ip3,1)+p(ip4,1))**2- # (p(ip3,2)+p(ip4,2))**2-(p(ip3,3)+p(ip4,3))**2 call lushow(ip3,ip4,sqrt(max(0.,pm34s))) endif * * do fragmentation and decays. possibly except tau decay. * if(itau.eq.0) then if(iabs(k(i1,2)).eq.15) k(i1,1)= 11 if(iabs(k(i2,2)).eq.15) k(i2,1)= 11 if(iabs(k(i3,2)).eq.15) k(i3,1)= 11 if(iabs(k(i4,2)).eq.15) k(i4,1)= 11 endif call luexec if(itau.eq.0) then if(iabs(k(i1,2)).eq.15) k(i1,1)= 1 if(iabs(k(i2,2)).eq.15) k(i2,1)= 1 if(iabs(k(i3,2)).eq.15) k(i3,1)= 1 if(iabs(k(i4,2)).eq.15) k(i4,1)= 1 endif * end * *-----wtoregion----------------------------------------------------- * subroutine wtoregion(n,x,j,a,b) implicit real*8(a-h,o-z) * dimension x(n) * if(j.gt.1) then do i=1,j-1 if(x(i).lt.0.d0.or.x(i).gt.1.d0) then stop endif enddo endif * a= 0.d0 b= 1.d0 * return end * *-------------------------------------------------------------------------- * subroutine wtocorrqcd(scal,als,bqm2,cqm2,rbm2,rcm2,rsm2) implicit real*8(a-h,o-z) * common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi * data rz3/1.20205690315959428540d0/,nf/5/ * scal2= scal*scal bqm= sqrt(bqm2) cqm= sqrt(cqm2) * *-----computes the running charm mass *-----computes the running bottom mass * scal1g= 1.d0 als1g= wtoralphas(wm,scal1g,als,nf)/pi alsc4= wtoralphas(wm,cqm,als,nf)/pi alsc42= alsc4*alsc4 als1g2= als1g*als1g alsb= wtoralphas(wm,bqm,als,nf)/pi alsb2= alsb*alsb alsz= wtoralphas(wm,scal,als,nf)/pi alsz2= alsz*alsz * *-----first the running of the s-quark mass * up to c/b-threshold * fn= 3.d0 b03= (11.d0-2.d0/3.d0*fn)/4.d0 b13= (102.d0-38.d0/3.d0*fn)/16.d0 b23= (2857.d0/2.d0-5033.d0/18.d0*fn+325.d0/54.d0* # fn*fn)/64.d0 g03= 1.d0 g13= (202.d0/3.d0-20.d0/9.d0*fn)/16.d0 g23= (1249.d0-(2216.d0/27.d0+160.d0/3.d0*rz3)*fn- # 140.d0/81.d0*fn*fn)/64.d0 rex3= g03/b03 b03s= b03*b03 b03c= b03s*b03 b13s= b13*b13 sm1g= 0.189d0 alsdf= alsc4-als1g cfm13= g13/b03-b13*g03/b03s cfm23= g23/b03-b13*g13/b03s-b23*g03/b03s+ # b13s*g03/b03c rsmc= sm1g*(alsc4/als1g)**rex3*(1.d0+cfm13* # alsdf+0.5d0*cfm13*cfm13*alsdf*alsdf+0.5d0*cfm23* # (alsc42-als1g2)) * fn= 4.d0 b04= (11.d0-2.d0/3.d0*fn)/4.d0 b14= (102.d0-38.d0/3.d0*fn)/16.d0 b24= (2857.d0/2.d0-5033.d0/18.d0*fn+325.d0/54.d0* # fn*fn)/64.d0 g04= 1.d0 g14= (202.d0/3.d0-20.d0/9.d0*fn)/16.d0 g24= (1249.d0-(2216.d0/27.d0+160.d0/3.d0*rz3)*fn- # 140.d0/81.d0*fn*fn)/64.d0 rex4= g04/b04 b04s= b04*b04 b04c= b04s*b04 b14s= b14*b14 alsdf= alsb-alsc4 cfm14= g14/b04-b14*g04/b04s cfm24= g24/b04-b14*g14/b04s-b24*g04/b04s+ # b14s*g04/b04c rsmb= rsmc*(alsb/alsc4)**rex4*(1.d0+cfm14* # alsdf+0.5d0*cfm14*cfm14*alsdf*alsdf+0.5d0*cfm24* # (alsb2-alsc42)) * *-----c quark mass at c-threshold * zero= 0.d0 x1= 0.5d0 x2= 2.0d0 xacc= 1.d-12 fn= 4.d0 cmm4= wtorrunm(x1,x2,xacc,cqm,alsc4,rsmc,zero,fn) * *-----c quark mass at b-threshold * cmb= cmm4*(alsb/alsc4)**rex4*(1.d0+cfm14* # (alsb-alsc4)+0.5d0*cfm14*cfm14* # (alsb-alsc4)**2+0.5d0*cfm24*(alsb2-alsc42)) * *-----running c mass * fn= 5.d0 b05= (11.d0-2.d0/3.d0*fn)/4.d0 b15= (102.d0-38.d0/3.d0*fn)/16.d0 b25= (2857.d0/2.d0-5033.d0/18.d0*fn+325.d0/54.d0* # fn*fn)/64.d0 g05= 1.d0 g15= (202.d0/3.d0-20.d0/9.d0*fn)/16.d0 g25= (1249.d0-(2216.d0/27.d0+160.d0/3.d0*rz3)*fn- # 140.d0/81.d0*fn*fn)/64.d0 rex5= g05/b05 b05s= b05*b05 b05c= b05s*b05 b15s= b15*b15 cfm15= g15/b05-b15*g05/b05s cfm25= g25/b05-b15*g15/b05s-b25*g05/b05s+b15s*g05/b05c rcm= cmb*(alsz/alsb)**rex5*(1.d0+cfm15* # (alsz-alsb)+0.5d0*cfm15*cfm15* # (alsz-alsb)**2+0.5d0*cfm25*(alsz2-alsb2)) rcm2= rcm*rcm * *-----running s mass * rsm= rsmb*(alsz/alsb)**rex5*(1.d0+cfm15* # (alsz-alsb)+0.5d0*cfm15*cfm15* # (alsz-alsb)**2+0.5d0*cfm25*(alsz2-alsb2)) rsm2= rsm*rsm * *-----b quark mass at b-threshold * x1= 0.5d0 x2= 6.0d0 xacc= 1.d-12 fn= 5.d0 bmm5= wtorrunm(x1,x2,xacc,bqm,alsb,cmb,rsmb,fn) * *-----running b mass * alsdf= alsz-alsb rbm= bmm5*(alsz/alsb)**rex5*(1.d0+cfm15* # alsdf+0.5d0*cfm15*cfm15*alsdf*alsdf+0.5d0*cfm25* # (alsz2-alsb2)) rbm2= rbm*rbm * return end * *-----wtoralphas-------------------------------------------------------------- *-----computes alpha_s(rs) for nf flavors given als= alpha_s(wm) * real*8 function wtoralphas(rs0,rs,als,nf) implicit real*8(a-h,o-z) * common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi common/wtfmass/em,rmm,tm,rnm,uqm,dqm,cqm,sqm,bqm,tqm,dmy * dimension b0(5),b1(5),b2(5) * *-----limits for lambda_5 are 1 mev < lambda_5 < 10 gev * if(als.eq.0.d0) then wtoralphas= 0.d0 else x1= 0.001d0 x2= 10.0d0 xacc= 1.d-12 qcdl= wtoqcdlam(nf,als,rs0,x1,x2,xacc) pqcdl5= qcdl do i=1,5 b0(i)= (11.d0-2.d0/3.d0*i)/4.d0 b1(i)= (102.d0-38.d0/3.d0*i)/16.d0 b2(i)= 0.5d0*(2857.d0-i*(5033.d0/9.d0- # 325.d0/27.d0*i))/64.d0 enddo * if(rs.lt.bqm) then rat= bqm/qcdl rl= 2.d0*log(rat) rll= log(rl) rb= log(b0(5)/b0(4)) fac= b1(5)/b0(5)-b1(4)/b0(4) facp= b2(5)/b0(5)-b2(4)/b0(4) fac2= b1(5)*b1(5)/b0(5)/b0(5)-b1(4)*b1(4)/b0(4)/b0(4) rhs= (b0(5)-b0(4))*rl+fac*rll-b1(4)/b0(4)*rb+ # b1(5)/b0(5)/b0(5)*fac*rll/rl+1.d0/b0(5)/rl*( # fac2-facp-7.d0/72.d0) rhs= rhs/b0(4) ratl2= exp(rhs) qcdl= qcdl/sqrt(ratl2) pqcdl4= qcdl nfe= nf-1 if(rs.lt.cqm) then rat= cqm/qcdl rl= 2.d0*log(rat) rll= log(rl) rb= log(b0(4)/b0(3)) fac= b1(4)/b0(4)-b1(3)/b0(3) facp= b2(4)/b0(4)-b2(3)/b0(3) fac2= b1(4)*b1(4)/b0(4)/b0(4)-b1(3)*b1(3)/b0(3)/b0(3) rhs= (b0(4)-b0(3))*rl+fac*rll-b1(3)/b0(3)*rb+ # b1(4)/b0(4)/b0(4)*fac*rll/rl+1.d0/b0(4)/rl*( # fac2-facp-7.d0/72.d0) rhs= rhs/b0(3) ratl2= exp(rhs) qcdl= qcdl/sqrt(ratl2) pqcdl3= qcdl nfe= nf-2 endif else nfe= nf endif * qcdb0= 11.d0-2.d0/3.d0*nfe qcdb1= 102.d0-38.d0/3.d0*nfe qcdb2= 0.5d0*(2857.d0-5033.d0/9.d0*nfe+ # 325.d0/27.d0*nfe*nfe) qcda= 2.d0*log(rs/qcdl) * wtoralphas= 4.d0*pi/qcdb0/qcda*(1.d0-qcdb1/qcdb0**2/qcda* # log(qcda)+(qcdb1/qcdb0**2/qcda)**2*((log(qcda)- # 0.5d0)**2+qcdb2*qcdb0/qcdb1**2-5.d0/4.d0)) * endif * return end * *-------------------------------------------------------------------- * real*8 function wtorals(qcdl,rs,nf) implicit real*8(a-h,o-z) * common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi common/wtfmass/em,rmm,tm,rnm,uqm,dqm,cqm,sqm,bqm,tqm,dmy * dimension b0(5),b1(5),b2(5) * pqcdl5= qcdl do i=1,5 b0(i)= (11.d0-2.d0/3.d0*i)/4.d0 b1(i)= (102.d0-38.d0/3.d0*i)/16.d0 b2(i)= 0.5d0*(2857.d0-i*(5033.d0/9.d0- # 325.d0/27.d0*i))/64.d0 enddo * if(rs.lt.bqm) then rat= bqm/qcdl rl= 2.d0*log(rat) rll= log(rl) rb= log(b0(5)/b0(4)) fac= b1(5)/b0(5)-b1(4)/b0(4) facp= b2(5)/b0(5)-b2(4)/b0(4) fac2= b1(5)*b1(5)/b0(5)/b0(5)-b1(4)*b1(4)/b0(4)/b0(4) rhs= (b0(5)-b0(4))*rl+fac*rll-b1(4)/b0(4)*rb+ # b1(5)/b0(5)/b0(5)*fac*rll/rl+1.d0/b0(5)/rl*( # fac2-facp-7.d0/72.d0) rhs= rhs/b0(4) ratl2= exp(rhs) qcdl= qcdl/sqrt(ratl2) pqcdl4= qcdl nfe= nf-1 if(rs.lt.cqm) then rat= cqm/qcdl rl= 2.d0*log(rat) rll= log(rl) rb= log(b0(4)/b0(3)) fac= b1(4)/b0(4)-b1(3)/b0(3) facp= b2(4)/b0(4)-b2(3)/b0(3) fac2= b1(4)*b1(4)/b0(4)/b0(4)-b1(3)*b1(3)/b0(3)/b0(3) rhs= (b0(4)-b0(3))*rl+fac*rll-b1(3)/b0(3)*rb+ # b1(4)/b0(4)/b0(4)*fac*rll/rl+1.d0/b0(4)/rl*( # fac2-facp-7.d0/72.d0) rhs= rhs/b0(3) ratl2= exp(rhs) qcdl= qcdl/sqrt(ratl2) pqcdl3= qcdl nfe= nf-2 endif else nfe= nf endif * qcdb0= 11.d0-2.d0/3.d0*nfe qcdb1= 102.d0-38.d0/3.d0*nfe qcdb2= 0.5d0*(2857.d0-5033.d0/9.d0*nfe+ # 325.d0/27.d0*nfe*nfe) qcda= 2.d0*log(rs/qcdl) * wtorals= 4.d0*pi/qcdb0/qcda*(1.d0-qcdb1/qcdb0**2/qcda* # log(qcda)+(qcdb1/qcdb0**2/qcda)**2*((log(qcda)- # 0.5d0)**2+qcdb2*qcdb0/qcdb1**2-5.d0/4.d0)) * return end * *-----wtorrunm-------------------------------------------------------- * computes the running quark mass at the pole mass * real*8 function wtorrunm(x1,x2,xacc,qm,als,rm1,rm2,fn) implicit real*8(a-h,o-z) * parameter (jmax=50) * fmid= wtoqcdmass(qm,als,rm1,rm2,fn,x2) f= wtoqcdmass(qm,als,rm1,rm2,fn,x1) if (f*fmid.ge.0.d0) then print*,'root must be bracketed for bisection' print 1,qm 1 format(/' error detected by rrunm ',/ # ' current value of quark mass = ',e20.5) stop endif if (f.lt.0.d0) then wtorrunm= x1 dx= x2-x1 else wtorrunm= x2 dx= x1-x2 endif do j=1,jmax dx= dx*0.5d0 xmid= wtorrunm+dx fmid= wtoqcdmass(qm,als,rm1,rm2,fn,xmid) if (fmid.le.0.d0) wtorrunm= xmid if(abs(dx).lt.xacc.or.fmid.eq.0.d0) return enddo pause 'too many bisections' end * *-----wtoqcdmass------------------------------------------------------ * real*8 function wtoqcdmass(qm,als,rm1,rm2,fn,x) implicit real*8(a-h,o-z) * common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi * rz2= pis/6.d0 rz3= 1.20205690315959428540d0 ra4= 0.5174790617d0 rz5= 1.03692775514336992633d0 rln= 2.d0*log(qm/x) rlns= rln*rln r1= rm1/x r2= rm2/x delta0= 3.d0/4.d0*rz2-3.d0/8.d0 delta1= r1*(pis/8.d0+r1*(-0.597d0+0.230d0*r1)) delta2= r2*(pis/8.d0+r2*(-0.597d0+0.230d0*r2)) * rhs= 1.d0+als*(4.d0/3.d0+rln+als*(3817.d0/288.d0-8.d0/3.d0+ # 2.d0/3.d0*(2.d0+log(2.d0))*rz2-rz3/6.d0-fn/3.d0*(rz2+ # 71.d0/48.d0)+4.d0/3.d0*(delta0+delta1+delta2)+(173.d0/ # 24.d0-13.d0/36.d0*fn)*rln+(15.d0/8.d0-fn/12.d0)*rlns)) wtoqcdmass= qm-x*rhs * return end * *-----wtoqcdlam----------------------------------------------------------------- *-----computes lamba^nf_msbar from alpha_s(rs0) * real*8 function wtoqcdlam(nf,als,rs,x1,x2,xacc) implicit real*8(a-h,o-z) * parameter (jmax=50,nout=21) * fmid= wtoqcdscale(nf,als,rs,x2) f= wtoqcdscale(nf,als,rs,x1) if (f*fmid.ge.0.d0) then write(nout,fmt=*) 'root must be bracketed for bisection' print*,'root must be bracketed for bisection' write(nout,fmt=1) als print 1,als 1 format(/' error detected by wtoqcdlam ',/ # ' current value of alpha_s = ',e20.5) stop endif if (f.lt.0.d0) then wtoqcdlam= x1 dx= x2-x1 else wtoqcdlam= x2 dx= x1-x2 endif do j=1,jmax dx= dx*0.5d0 xmid= wtoqcdlam+dx fmid= wtoqcdscale(nf,als,rs,xmid) if (fmid.le.0.d0) wtoqcdlam= xmid if(abs(dx).lt.xacc.or.fmid.eq.0.d0) return enddo pause 'too many bisections' end * *-----wtoqcdscale------------------------------------------------------------- *-----computes lamba^nf_msbar from alpha_s(rs0) * real*8 function wtoqcdscale(nf,als,rs,x) implicit real*8(a-h,o-z) * common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi * qcdb0= 11.d0-2.d0/3.d0*nf qcdb1= 102.d0-38.d0/3.d0*nf qcdb2= 0.5d0*(2857.d0-5033.d0/9.d0*nf+325.d0/27.d0*nf*nf) qcda= 2.d0*log(rs/x) wtoqcdscale= als-(4.d0*pi/qcdb0/qcda*(1.d0-qcdb1/qcdb0**2/qcda* # log(qcda)+(qcdb1/qcdb0**2/qcda)**2*((log(qcda)- # 0.5d0)**2+qcdb2*qcdb0/qcdb1**2-5.d0/4.d0))) * return end * subroutine wtohadr5(e,der,eder) ****************************************************************** * * * subroutine for the evaluation of the light hadron * * contributions to delta_r and delta_g * * using fits to the * * qed vacuum polarization from e^+ e^- data * * * * f. jegerlehner, paul scherrer institute, ch-5232 villigen * * * * reference: f. jegerlehner, z. phys. c32 (1986) 195 * * h. burkhardt et al., z. phys. c42 (1989) 497 * * s. eidelman, f. jegerlehner, z. phys. c (1995) * * * ****************************************************************** * version: 24/02/1995 * * notation: e energy ( momentum transfer ): e>0 timelike , e<0 spacelike * st2 is sin^2(theta); st2= 0.2322 is the reference value * the routine returns the hadronic contribution of 5 flavors (u,d,s,c,b) * to der= delta_r with hadronic error errder * and deg= delta_g with hadronic error errdeg * the effective value of the fine structure constant alphaqed at energy * e is alphaqed(e)= alphaqed(0)/(1-delta_r) ,similarly for the su(2) * coupling alphasu2(e)= alphasu2(0)/(1-delta_g), where delta_r(g) is the * sum of leptonic, hadronic contributions (top to be added). * * this program does not yet know how to compute delta r and delta g for * energies in the ranges |e|>1tev and 2m_pi < e < 40(13) gev !!!!!!!!! * implicit real*8(a-h,o-z) parameter(nf=9,ns=4) * dimension res(ns) dimension ae(nf,ns),eu(nf),eo(nf),rm1(nf) dimension c1(nf,ns),c2(nf,ns),c3(nf,ns),c4(nf,ns),rl1(nf,ns) * do i=1,nf eu(i)= 0.d0 eo(i)= 0.d0 rm1(i)= 0.d0 do j=1,ns ae(i,j)= 0.d0 c1(i,j)= 0.d0 c2(i,j)= 0.d0 c3(i,j)= 0.d0 c4(i,j)= 0.d0 enddo enddo * * #1# delta_r * fit parameters spacelike -1000 to -200 gev eu(1) = -1.d3 eo(1) = -2.d2 rm1(1) = -1.d3 c1(1,1)= 4.2069394d-02 c2(1,1)= 2.9253566d-03 c3(1,1)= -6.7782454d-04 c4(1,1)= 9.3214130d-06 * fit parameters spacelike -200 to -20 gev eu(2) = -2.d2 eo(2) = -2.d1 rm1(2) = -1.d2 c1(2,1)= 2.8526291d-02 c2(2,1)= 2.9520725d-03 c3(2,1)= -2.7906310d-03 c4(2,1)= 6.4174528d-05 * fit parameters spacelike -20 to -2 gev eu(3) = -2.d1 eo(3) = -2.d0 rm1(3) = -2.d1 rl1(3,1)= 9.3055d-03 c1(3,1)= 2.8668314d-03 c2(3,1)= 0.3514608d0 c3(3,1)= 0.5496359d0 c4(3,1)= 1.9892334d-04 ae(3,1)= 3.d0 * fit parameters spacelike -2 to 0.25 gev eu(4) = -2.d0 eo(4) = 0.25d0 rm1(4) = -2.d0 rl1(4,1)= 9.3055d-03 c1(4,1)= 2.2694240d-03 c2(4,1)= 8.073429d0 c3(4,1)= 0.1636393d0 c4(4,1)= -3.3545541d-05 ae(4,1)= 2.d0 * fit parameters timelike 0.25 to 2 gev eu(5) = 0.25d0 eo(5) = 2.d0 * fit parameters timelike 2 to 40 gev eu(6) = 2.d0 eo(6) = 4.d1 * fit parameters timelike 40 to 80 gev eu(7) = 4.d1 eo(7) = 8.d1 rm1(7) = 8.d1 c1(7,1)= 2.7266588d-02 c2(7,1)= 2.9285045d-03 c3(7,1)= -4.7720564d-03 c4(7,1)= 7.7295507d-04 * fit parameters timelike 80 to 250 gev eu(8) = 8.d1 eo(8) = 2.5d2 rm1(8) = 91.18880d0 c1(8,1)= 2.8039809d-02 c2(8,1)= 2.9373798d-03 c3(8,1)= -2.8432352d-03 c4(8,1)= -5.2537734d-04 * fit parameters timelike 250 to 1000 gev eu(9) = 2.5d2 eo(9) = 1.d3 rm1(9) = 1.d3 c1(9,1)= 4.2092260d-02 c2(9,1)= 2.9233438d-03 c3(9,1)= -3.2966913d-04 c4(9,1)= 3.4324117d-07 * #2# delta_g * fit parameters spacelike -1000 to -200 gev c1(1,2)= 8.6415343d-02 c2(1,2)= 6.0127582d-03 c3(1,2)= -6.7379221d-04 c4(1,2)= 9.0877611d-06 * fit parameters spacelike -200 to -20 gev c1(2,2)= 5.8580618d-02 c2(2,2)= 6.0678599d-03 c3(2,2)= -2.4153464d-03 c4(2,2)= 6.1934326d-05 * fit parameters spacelike -20 to -2 gev rl1(3,2)= 1.9954d-02 c1(3,2)= 5.7231588d-03 c2(3,2)= 0.3588257d0 c3(3,2)= 0.5532265d0 c4(3,2)= 6.0730567d-04 ae(3,2)= 3.d0 * fit parameters spacelike -2 to 0.25 gev rl1(4,2)= 1.9954d-02 c1(4,2)= 4.8065037d-03 c2(4,2)= 8.255167d0 c3(4,2)= 0.1599882d0 c4(4,2)= -1.8624817d-04 ae(3,2)= 2.d0 * fit parameters timelike 40 to 80 gev c1(7,2)= 5.5985276d-02 c2(7,2)= 6.0203830d-03 c3(7,2)= -5.0066952d-03 c4(7,2)= 7.1363564d-04 * fit parameter timelike 80 to 250 gev c1(8,2)= 5.7575710d-02 c2(8,2)= 6.0372148d-03 c3(8,2)= -3.4556778d-03 c4(8,2)= -4.9574347d-04 * fit parameters timelike 250 to 1000 gev c1(9,2)= 8.6462371d-02 c2(9,2)= 6.0088057d-03 c3(9,2)= -3.3235471d-04 c4(9,2)= 5.9021050d-07 * #3# error delta_r * fit parameters spacelike -1000 to -200 gev c1(1,3)= 6.3289929d-04 c2(1,3)= 3.3592437d-06 c3(1,3)= 0.d0 c4(1,3)= 0.d0 * fit parameters spacelike -200 to -20 gev c1(2,3)= 6.2759849d-04 c2(2,3)= -1.0816625d-06 c3(2,3)= 5.050189d0 c4(2,3)= -9.6505374d-02 ae(2,3)= 1.d0 * fit parameters spacelike -20 to -2 gev rl1(3,3)= 2.0243d-04 c1(3,3)= 1.0147886d-04 c2(3,3)= 1.819327d0 c3(3,3)= -0.1174904d0 c4(3,3)= -1.2404939d-04 ae(3,3)= 3.d0 * fit parameters spacelike -2 to 0.25 gev rl1(4,3)= 2.0243d-04 c1(4,3)= -7.1368617d-05 c2(4,3)= 9.980347d-04 c3(4,3)= 1.669151d0 c4(4,3)= 3.5645600d-05 ae(4,3)= 2.d0 * fit parameters timelike 40 to 80 gev c1(7,3)= 6.4947648d-04 c2(7,3)= 4.9386853d-07 c3(7,3)= -55.22332d0 c4(7,3)= 26.13011d0 * fit parameters timelike 80 to 250 gev c1(8,3)= 6.4265809d-04 c2(8,3)= -2.8453374d-07 c3(8,3)= -23.38172d0 c4(8,3)= -6.251794d0 * fit parameter timelike 250 to 1000 gev c1(9,3)= 6.3369947d-04 c2(9,3)= -2.0898329d-07 c3(9,3)= 0.d0 c4(9,3)= 0.d0 * #4# error delta_g * fit parameters spacelike -1000 to -200 gev c1(1,4)= 1.2999176d-03 c2(1,4)= 7.4505529d-06 c3(1,4)= 0.d0 c4(1,4)= 0.d0 * fit parameters spacelike -200 to -20 gev c1(2,4)= 1.2883141d-03 c2(2,4)= -1.3790827d-06 c3(2,4)= 8.056159d0 c4(2,4)= -0.1536313d0 ae(2,4)= 1.d0 * fit parameters spacelike -20 to -2 gev rl1(3,4)= 4.3408d-04 c1(3,4)= 2.0489733d-04 c2(3,4)= 2.065011d0 c3(3,4)= -0.6172962d0 c4(3,4)= -2.5603661d-04 ae(3,4)= 3.d0 * fit parameters spacelike -2 to 0.25 gev rl1(4,4)= 4.3408d-04 c1(4,4)= -1.5095409d-04 c2(4,4)= 9.9847501d-04 c3(4,4)= 1.636659d0 c4(4,4)= 7.5892596d-05 ae(4,4)= 2.d0 * fit parameters timelike 40 to 80 gev c1(7,4)= 1.3335156d-03 c2(7,4)= 2.2939612d-07 c3(7,4)= -246.4966d0 c4(7,4)= 114.9956d0 * fit parameters timelike 80 to 250 gev c1(8,4)= 1.3196438d-03 c2(8,4)= 2.8937683d-09 c3(8,4)= 5449.778d0 c4(8,4)= 930.3875d0 * fit parameters timelike 250 to 1000 gev c1(9,4)= 1.3016918d-03 c2(9,4)= -3.6027674d-07 c3(9,4)= 0.d0 c4(9,4)= 0.d0 * se= 654.d0/643.d0 ! rescaling error to published version 1995 s= e*e der= 0.d0 deg= 0.d0 eder= 0.d0 edeg= 0.d0 if((e.gt.1.e3).or.(e.lt.-1.e3)) goto 100 if((e.lt.4.e1).and.(e.gt.0.25d0)) goto 100 i= 1 do while (e.ge.eo(i)) i= i+1 enddo if(e.eq.1.e3) i= 9 if(e.eq.0.d0 ) goto 100 s0= sign(1.d0,rm1(i))*rm1(i)**2 s = sign(1.d0,e)*e*e x1= s0/s xi= 1.d0/x1 x2= x1*x1 if(ae(i,1).le.0.d0) then do j= 1,4 xlar= xi+ae(i,j)*exp(-xi) xlog= log(xlar) res(j)= c1(i,j)+c2(i,j)*(xlog+c3(i,j)*(x1-1.d0)+ # c4(i,j)*(x2-1.0)) enddo else if (ae(i,1).eq.2.d0) then hx= xi*xi do j= 1,2 fx= 1.d0-c2(i,j)*s gx= c3(i,j)*s/(c3(i,j)-s) xx= log(abs(fx))+c2(i,j)*gx res(j)= c1(i,j)*xx-rl1(i,j)*gx+c4(i,j)*hx enddo do j= 3,4 u= abs(s) gx= -c3(i,j)*u/(c3(i,j)+u) xx= xi**3/(sqrt(abs(xi))**5+c2(i,j)) res(j)= c1(i,j)*xx-rl1(i,j)*gx+c4(i,j)*hx enddo else if (ae(i,1).eq.3.0) then hx= xi do j= 1,4 fx= 1.d0-c2(i,j)*s gx= c3(i,j)*s/(c3(i,j)-s) xx= log(abs(fx))+c2(i,j)*gx res(j)= c1(i,j)*xx-rl1(i,j)*gx+c4(i,j)*hx enddo endif der= res(1) eder= res(3)*se 100 return end * *------------------------------------------------------------- * subroutine wtopself(p2x,pggf) implicit real*8(a-h,o-z) * common/wtparam/eps,ddelta common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi common/wtfmass/em,rmm,tm,rnm,uqm,dqm,cqm,sqm,bqm,tqm,dmy * em2= em*em amm2= rmm*rmm tm2= tm*tm tqm2= tqm*tqm * call wtorbff(p2x,em2,em2,b0ee,b1ee,b21ee) call wtorbff(p2x,amm2,amm2,b0mm,b1mm,b21mm) call wtorbff(p2x,tm2,tm2,b0tau,b1tau,b21tau) call wtorbff(p2x,tqm2,tqm2,b0tt,b1tt,b21tt) call wtorbff0(em2,em2,b0ee0,b1ee0,b21ee0) call wtorbff0(amm2,amm2,b0mm0,b1mm0,b21mm0) call wtorbff0(tm2,tm2,b0tau0,b1tau0,b21tau0) call wtorbff0(tqm2,tqm2,b0tt0,b1tt0,b21tt0) * bfe= 2.d0*b21ee-b0ee bfm= 2.d0*b21mm-b0mm bftau= 2.d0*b21tau-b0tau bft= 2.d0*b21tt-b0tt bfe0= 2.d0*b21ee0-b0ee0 bfm0= 2.d0*b21mm0-b0mm0 bftau0= 2.d0*b21tau0-b0tau0 bft0= 2.d0*b21tt0-b0tt0 * sllc= 1.d0+3.d0/4.d0/alphai/pi pggf= 4.d0*(bfe+bfm+bftau)*sllc+16.d0/3.d0*bft- # 4.d0*(bfe0+bfm0+bftau0)*sllc-16.d0/3.d0*bft0 * return end * *------------------------------------------------------------------ * subroutine wtopselfnp(op2x,opggnp) implicit real*8(a-h,o-z) * common/wtparam/eps,ddelta common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi common/wtfmass/em,rmm,tm,rnm,uqm,dqm,cqm,sqm,bqm,tqm,dmy * em2= em*em amm2= rmm*rmm tm2= tm*tm tqm2= tqm*tqm sllc= 1.d0+3.d0/4.d0/alphai/pi * p2x= op2x p2z= -zm*zm qzm2= -p2z call wtorbff(p2z,em2,em2,b0eez,b1eez,b21eez) call wtorbff(p2z,amm2,amm2,b0mmz,b1mmz,b21mmz) call wtorbff(p2z,tm2,tm2,b0tauz,b1tauz,b21tauz) call wtorbff(p2z,tqm2,tqm2,b0ttz,b1ttz,b21ttz) call wtorbff0(em2,em2,b0ee0,b1ee0,b21ee0) call wtorbff0(amm2,amm2,b0mm0,b1mm0,b21mm0) call wtorbff0(tm2,tm2,b0tau0,b1tau0,b21tau0) call wtorbff0(tqm2,tqm2,b0tt0,b1tt0,b21tt0) bfez= 2.d0*b21eez-b0eez bfmz= 2.d0*b21mmz-b0mmz bftauz= 2.d0*b21tauz-b0tauz bftz= 2.d0*b21ttz-b0ttz bfe0= 2.d0*b21ee0-b0ee0 bfm0= 2.d0*b21mm0-b0mm0 bftau0= 2.d0*b21tau0-b0tau0 bft0= 2.d0*b21tt0-b0tt0 pggl= 4.d0*(bfez+bfmz+bftauz)*sllc+16.d0/3.d0*bftz pggl0= 4.d0*(bfe0+bfm0+bftau0)*sllc+16.d0/3.d0*bft0 b= 0.00299d0 c= 1.d0 an= -b*log(1.d0+c*qzm2)-0.25d0/alphai/pi*(pggl-pggl0)+ # 1.d0-128.89d0/alphai da= an-0.00165d0 ap2x= abs(p2x) px= sqrt(ap2x) if(px.lt.0.3d0) then a= 0.d0 b= 0.00835d0 c= 1.d0 else if(px.lt.3.d0) then a= 0.d0 b= 0.00238d0 c= 3.927d0 else if(px.lt.100.d0) then a= 0.00165d0 b= 0.00299d0 c= 1.d0 else if(px.gt.100.d0) then a= 0.00221d0 b= 0.00293d0 c= 1.d0 endif * a= a+da opggnp= 4.d0*pi*alphai*(a+b*log(1.d0+c*ap2x)) * return end * *------------------------------------------------------------------- * subroutine wtorbff0(rm12,rm22,b0,b1,b21) implicit real*8(a-h,o-z) * common/wtparam/eps,ddelta common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi * dimension arg(2),cln(2),fr(3) * if(rm12.eq.0.d0) then arm12= 1.d-20 else arm12= rm12 endif if(rm22.eq.0.d0) then arm22= 1.d-20 else arm22= rm22 endif if (arm12.eq.arm22) then fact= ddelta-log(arm12) b0= fact b1= -0.5d0*fact b21= fact/3.d0 return else n= 3 yr= arm12/(arm12-arm22) omyr= -arm22/(arm12-arm22) yi= -eps/(arm12-arm22) call wtorcg(n,yr,yi,omyr,fr) arg(1)= omyr arg(2)= -yi call wtoclnx(arg,cln) f1r= fr(1)-cln(1) f2r= 2.d0*fr(2)-cln(1) f3r= 3.d0*fr(3)-cln(1) fact= ddelta-log(arm22) b0= fact-f1r b1= 0.5d0*(-fact+f2r) b21= 1.d0/3.d0*(fact-f3r) return endif end * *---------------------------------------------------------- * subroutine wtorbff(p2r,rm12,rm22,rb0,rb1,rb21) implicit real*8(a-h,o-z) * common/wtparam/eps,ddelta common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi * dimension cmp2(2) dimension clnmp2(2) dimension gfpr(3),gfmr(3) * if(rm12.eq.0.d0.and.rm22.eq.0.d0) then cmp2(1)= p2r cmp2(2)= -eps call wtoclnx(cmp2,clnmp2) rb0= ddelta+2.d0-clnmp2(1) rb1= -0.5d0*rb0 rb21= 1.d0/18.d0+1.d0/3.d0*rb0 else cmp2(1)= -p2r cmp2(2)= -eps call wtoroots(p2r,rm12,rm22,rpr,rpi,rmr,rmi,omrpr,omrmr) call wtoclnx(cmp2,clnmp2) n= 3 call wtorcg(n,rpr,rpi,omrpr,gfpr) call wtorcg(n,rmr,rmi,omrmr,gfmr) auxdel= ddelta-clnmp2(1) rb0= auxdel-gfpr(1)-gfmr(1) rb1= -auxdel/2.d0+gfpr(2)+gfmr(2) rb21= auxdel/3.d0-gfpr(3)-gfmr(3) endif * return end * *---------------------------------------------------------- * subroutine wtoclnx(arg,res) implicit real*8(a-h,o-z) * common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi * dimension arg(2),aarg(2),res(2) * do i= 1,2 aarg(i)= abs(arg(i)) enddo azm2= (arg(1))**2+(arg(2))**2 azm= sqrt(azm2) res(1)= log(azm) if (arg(1).eq.0.d0) then if (arg(2).gt.0.d0) then teta= pi/2.d0 else teta= -pi/2.d0 endif res(2)= teta return else if (arg(2).eq.0.d0) then if (arg(1).gt.0.d0) then teta= 0.d0 else teta= pi endif res(2)= teta return else tnteta= aarg(2)/aarg(1) teta= atan(tnteta) sr= arg(1)/aarg(1) si= arg(2)/aarg(2) if (sr.gt.0.d0) then res(2)= si*teta else res(2)= si*(pi-teta) endif return endif end * *--------------------------------------------------------- * subroutine wtoclnomx(arg,omarg,res) implicit real*8(a-h,o-z) * dimension arg(2),omarg(2),res(2),ares(2),ct(10),sn(10) * zr= arg(1) zi= arg(2) azm2= zr*zr+zi*zi azm= sqrt(azm2) if (azm.lt.1.d-7) then ct(1)= zr/azm sn(1)= zi/azm do n=2,10 ct(n)= ct(1)*ct(n-1)-sn(1)*sn(n-1) sn(n)= sn(1)*ct(n-1)+ct(1)*sn(n-1) enddo ares(1)= ct(10)/10.d0 ares(2)= sn(10)/10.d0 do k=9,1,-1 ares(1)= ares(1)*azm+ct(k)/k ares(2)= ares(2)*azm+sn(k)/k enddo ares(1)= -ares(1)*azm ares(2)= -ares(2)*azm else call wtoclnx(omarg,ares) endif do i= 1,2 res(i)= ares(i) enddo * return end * *---------------------------------------------------------- * subroutine wtocqlnx(arg,res) implicit real*16(a-h,o-z) * common/wtqparam/qpi,qpis,qeps,qdelta * dimension arg(2),aarg(2),res(2) * do i= 1,2 aarg(i)= abs(arg(i)) enddo zm2= (arg(1))**2+(arg(2))**2 zm= sqrt(zm2) res(1)= log(zm) if (arg(1).eq.0.q0) then if (arg(2).gt.0.q0) then teta= qpi/2.q0 else teta= -qpi/2.q0 endif res(2)= teta return else if (arg(2).eq.0.q0) then if (arg(1).gt.0.q0) then teta= 0.q0 else teta= qpi endif res(2)= teta return else tnteta= aarg(2)/aarg(1) teta= atan(tnteta) sr= arg(1)/aarg(1) si= arg(2)/aarg(2) if (sr.gt.0.q0) then res(2)= si*teta else res(2)= si*(qpi-teta) endif return endif end * *------------------------------------------------------------- * subroutine wtocqlnomx(arg,omarg,res) implicit real*16(a-h,o-z) * common/wtqparam/qpi,qpis,qeps,qdelta * dimension arg(2),aarg(2),omarg(2),res(2),ares(2), # ct(10),sn(10) * zr= arg(1) zi= arg(2) zm2= zr*zr+zi*zi zm= sqrt(zm2) if (zm.lt.1.q-7) then ct(1)= zr/zm sn(1)= zi/zm do n=2,10 ct(n)= ct(1)*ct(n-1)-sn(1)*sn(n-1) sn(n)= sn(1)*ct(n-1)+ct(1)*sn(n-1) enddo ares(1)= ct(10)/10.q0 ares(2)= sn(10)/10.q0 do k=9,1,-1 ares(1)= ares(1)*zm+ct(k)/k ares(2)= ares(2)*zm+sn(k)/k enddo ares(1)= -ares(1)*zm ares(2)= -ares(2)*zm else yr= omarg(1) yi= omarg(2) ym2= yr*yr+yi*yi ym= sqrt(ym2) do i= 1,2 aarg(i)= abs(omarg(i)) enddo res(1)= log(ym) if(yr.eq.0.q0) then if(yi.gt.0.q0) then teta= qpi/2.q0 else teta= -qpi/2.q0 endif ares(2)= teta else if(yi.eq.0.q0) then if(yr.gt.0.q0) then teta= 0.q0 else teta= qpi endif ares(2)= teta else tnteta= aarg(2)/aarg(1) teta= atan(tnteta) sr= yr/aarg(1) si= yi/aarg(2) if(sr.gt.0.q0) then ares(2)= si*teta else ares(2)= si*(qpi-teta) endif endif endif do i= 1,2 res(i)= ares(i) enddo * return end * *-------------------------------------------------------------------- * subroutine wtorcg(n,zr,zi,omzr,gfr) implicit real*8(a-h,o-z) * common/wtparam/eps,ddelta common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi * dimension gfr(n) dimension a(4),ra(16),z(2),omz(2),oz(2),clnomz(2),clnoz(2) dimension ca(2,10),aux(2),zp(2),ct(16),sn(16) * data ra/-0.2d+0,0.666666666666667d-1, # -0.952380952380952d-2,-0.396825396825397d-3, # 0.317460317460317d-3,-0.132275132275132d-4, # -0.962000962000962d-5,0.105218855218855d-5, # 0.266488361726450d-6,-0.488745528428000d-7, # -0.675397500794000d-8,0.190720263471000d-8, # 0.153663007690000d-9,-0.679697905790000d-10, # -0.293683556000000d-11,0.228836696000000d-11/ * z(1)= zr z(2)= zi omz(1)= omzr omz(2)= -zi if (zr.eq.0.d0.and.zi.eq.0.d0) then do k=1,n gfr(k)= -1.d0/k/k enddo return else if (zr.eq.1.d0.and.zi.eq.0.d0) then a(1)= -1.d0 do j=2,4 a(j)= ((j-1.d0)*a(j-1)-1.d0/j)/j enddo do k=1,n gfr(k)= a(k) enddo return else zmod2= zr*zr+zi*zi zmod= sqrt(zmod2) call wtoclnx(omz,clnomz) * |z| < 4 if (zmod.lt.4.d0) then oz(1)= -z(1) oz(2)= -z(2) call wtoclnx(oz,clnoz) ca(1,1)= omz(1)*clnomz(1)-omz(2)*clnomz(2)+z(1)*clnoz(1)- # z(2)*clnoz(2)-1.d0 ca(2,1)= omz(1)*clnomz(2)+omz(2)*clnomz(1)+z(1)*clnoz(2)+ # z(2)*clnoz(1) if (n.eq.1) then gfr(1)= ca(1,1) return else do j= 2,n jm1= j-1 ca(1,j)= ((j-1.d0)*(z(1)*ca(1,jm1)-z(2)*ca(2,jm1))+ # omz(1)*clnomz(1)-omz(2)*clnomz(2)-1.d0/j)/j ca(2,j)= ((j-1.d0)*(z(1)*ca(2,jm1)+z(2)*ca(1,jm1))+ # omz(1)*clnomz(2)+omz(2)*clnomz(1))/j enddo do k=1,n gfr(k)= ca(1,k) enddo return endif * |z| > 4 else aux(1)= (-zr*omzr+zi*zi)/zmod2 aux(2)= zi/zmod2 call wtoclnx(aux,zp) zpm2= zp(1)*zp(1)+zp(2)*zp(2) zpm= sqrt(zpm2) ct(1)= zp(1)/zpm sn(1)= zp(2)/zpm do k=2,16 ct(k)= ct(1)*ct(k-1)-sn(1)*sn(k-1) sn(k)= sn(1)*ct(k-1)+ct(1)*sn(k-1) enddo ca(1,4)= ra(16)*ct(16)*zpm+ra(15)*ct(15) ca(2,4)= ra(16)*sn(16)*zpm+ra(15)*sn(15) do j=14,1,-1 ca(1,4)= ca(1,4)*zpm+ra(j)*ct(j) ca(2,4)= ca(2,4)*zpm+ra(j)*sn(j) enddo ca(1,4)= ca(1,4)*zpm ca(2,4)= ca(2,4)*zpm do j= 3,1,-1 jp1= j+1 ca(1,j)= ((ca(1,jp1)+1.d0/jp1)*z(1)+ca(2,jp1)*z(2))/zmod2 ca(2,j)= (ca(2,jp1)*z(1)-(ca(1,jp1)+1.d0/jp1)*z(2))/zmod2 enddo do k=1,n gfr(k)= (ca(1,k)+clnomz(1))/k enddo return endif endif end * * *-----rspence----------------------------------------------------- * gives the real part of li_2(z) in double precision. accuracy * is about 8 digits * real*8 function wtorspence(xr,xi) implicit real*8(a,b,d-h,o-z) implicit complex*16 (c) * data pi/3.141592653589793238462643d0/ * data b0/1.d0/,b1/-0.25d0/,b2/0.277777777777778d-1/, # b4/-0.277777777777778d-3/,b6/0.472411186696901d-5/ * c0= cmplx(0.d0) pis= pi*pi pis6= pis/6.d0 cx= cmplx(xr,xi) comx= 1.d0-cx if (xr.lt.0.d0) then cy= 1.d0-cx sign1= -1.d0 clnx= log(cx) clnomx= log(comx) cadd1= pis6-clnx*clnomx else cy= cx sign1= 1.d0 cadd1= c0 endif comy= 1.d0-cy ym2= conjg(cy)*cy ym= sqrt(ym2) if (ym.gt.1.d0) then cz= conjg(cy)/ym2 sign2= -1.d0 coy= -cy clnoy= log(coy) cadd2= -pis6-0.5d0*clnoy*clnoy else cz= cy sign2= 1.d0 cadd2= c0 endif comz= 1.d0-cz zr= real(cz) if (zr.gt.0.5d0) then ct= 1.d0-cz comt= 1.d0-ct sign3= -1.d0 clnz= log(cz) clnomz= log(comz) cadd3= pis6-clnz*clnomz else ct= cz comt= 1.d0-ct sign3= 1.d0 cadd3= c0 endif cpar= log(comt) cpar2= cpar*cpar * cres= -cpar*(b0-cpar*(b1-cpar*(b2+cpar2*(b4+b6*cpar2)))) * cli2= sign1*(sign2*(sign3*cres+cadd3)+cadd2)+cadd1 * wtorspence= real(cli2) * return end * *---------------------------------------------------------- * real*8 function wtorfsr(x,r) implicit real*8(a,b,d-h,o-z) implicit complex*16(c) * x2= x*x r2= r*r * a1r= 1.d0 a1i= 1.d0/r a2r= 1.d0/(1.d0+r2*x2) a2i= r*x/(1.d0+r2*x2) a3r= 1.d0/(1.d0+r2) a3i= r/(1.d0+r2) a4r= (r2*x+1.d0)/(1.d0+r2) a4i= r*(1.d0-x)/(1.d0+r2) * ca1= cmplx(a1r,a1i) ca2= cmplx(a2r,a2i) * cpar= -log(ca1)*log(ca2) * wtorfsr= 2.d0*(real(cpar)+wtorspence(a3r,a3i)- # wtorspence(a4r,a4i)) * return end * *-------------------------------------------------------------- * subroutine funct1(n,xc,fc) implicit real*8(a-h,o-z) character*1 omssm character*4 otype parameter(npos=512) * common/wtkr/ireg common/wtmssmo/omssm common/wtochannel/otype common/wthx0/xshmx0(npos) common/wtcx0/xscmx0(npos) * dimension xc(n) * if(otype.eq.'cc11'.or.otype.eq.'cc20') then fc= -wtoxsc(n,xc)/xscmx0(ireg) else if(otype.eq.'cc12') then fc= -wtoxssc(n,xc)/xshmx0(ireg) else if(otype.eq.'nc33') then fc= -wtoxsnh32(n,xc)/xshmx0(ireg) else if(otype.eq.'nc32') then fc= -wtoxsn32(n,xc)/xshmx0(ireg) else if(otype.eq.'nc50') then fc= -wtoxsnh49(n,xc)/xshmx0(ireg) else if(otype.eq.'nc25') then fc= -wtoxsnh24(n,xc)/xshmx0(ireg) else if(otype.eq.'nc21') then fc= -wtoxsnh19(n,xc)/xshmx0(ireg) else if(otype.eq.'nc19') then fc= -wtoxsn(n,xc)/xshmx0(ireg) else if(otype.eq.'nc24') then fc= -wtoxsn(n,xc)/xshmx0(ireg) else if(otype.eq.'nc26') then fc= -wtoxsnh26(n,xc)/xshmx0(ireg) else if(otype.eq.'nc48') then fc= -wtoxsn(n,xc)/xshmx0(ireg) else if(otype.eq.'nc64') then fc= -wtoxsn64(n,xc)/xshmx0(ireg) else if(otype.eq.'nc68') then if(omssm.eq.'n') then fc= -wtoxsnh64(n,xc)/xshmx0(ireg) else if(omssm.eq.'y') then fc= -wtoxsnsh64(n,xc)/xshmx0(ireg) endif endif * return end * *------------------------------------------------------ * integer*4 function iwtopos(n,x) implicit real*8(a-h,o-z) character*1,osm * common/wtsmod/osm common/wtspl/prx0(9) common/wtspln/prx0n(9) * dimension x(n),i(9),apx(9) * do j=1,9 if(osm.eq.'n') then apx(j)= prx0(j) else if(osm.eq.'g') then apx(j)= prx0n(j) endif enddo * do j=1,n if(x(j).lt.apx(j)) then i(j)= 1 else i(j)= 2 endif enddo * iwtopos= 2*(2*(2*(2*(2*(2*(2*(2*(i(1)-1)+i(2)-1)+ # i(3)-1)+i(4)-1)+i(5)-1)+i(6)-1)+i(7)-1)+i(8)-1)+i(9) * return end * *------------------------------------------------------------------------- * subroutine wtomxregion(ipos,bl,bu) implicit real*8(a-h,o-z) * common/wtspl/prx0(9) * dimension bl(9),bu(9) * data eps/1.d-6/ * do i1=1,2 do i2=1,2 do i3=1,2 do i4=1,2 do i5=1,2 do i6=1,2 do i7=1,2 do i8=1,2 do i9=1,2 ip= 2*(2*(2*(2*(2*(2*(2*(2*(i1-1)+i2-1)+ # i3-1)+i4-1)+i5-1)+i6-1)+i7-1)+i8-1)+i9 if(ip.eq.ipos) then if(i1.eq.1) then bl(1)= eps bu(1)= prx0(1) else bl(1)= prx0(1) bu(1)= 1.d0-eps endif if(i2.eq.1) then bl(2)= eps bu(2)= prx0(2) else bl(2)= prx0(2) bu(2)= 1.d0-eps endif if(i3.eq.1) then bl(3)= eps bu(3)= prx0(3) else bl(3)= prx0(3) bu(3)= 1.d0-eps endif if(i4.eq.1) then bl(4)= eps bu(4)= prx0(4) else bl(4)= prx0(4) bu(4)= 1.d0-eps endif if(i5.eq.1) then bl(5)= eps bu(5)= prx0(5) else bl(5)= prx0(5) bu(5)= 1.d0-eps endif if(i6.eq.1) then bl(6)= eps bu(6)= prx0(6) else bl(6)= prx0(6) bu(6)= 1.d0-eps endif if(i7.eq.1) then bl(7)= eps bu(7)= prx0(7) else bl(7)= prx0(7) bu(7)= 1.d0-eps endif if(i8.eq.1) then bl(8)= eps bu(8)= prx0(8) else bl(8)= prx0(8) bu(8)= 1.d0-eps endif if(i9.eq.1) then bl(9)= eps bu(9)= prx0(9) else bl(9)= prx0(9) bu(9)= 1.d0-eps endif endif enddo enddo enddo enddo enddo enddo enddo enddo enddo * return end * *---------------------------------------------------------------- * subroutine wtovolume(npos,nint,xsmx,oxsmx,gbl,gbu,vol,volt) implicit real*8(a-h,o-z) character*1 ord parameter(nposm=512) * common/wtspln/prx0n(9) * dimension vol(npos,nint),volt(nint) dimension avol(nposm),ir(nposm) dimension bl(nposm,9),bu(nposm,9) dimension xsmx(npos,nint),oxsmx(npos,nint),gbl(npos,9,nint), # gbu(npos,9,nint) * do i=1,npos do i1=1,2 do i2=1,2 do i3=1,2 do i4=1,2 do i5=1,2 do i6=1,2 do i7=1,2 do i8=1,2 do i9=1,2 ip= 2*(2*(2*(2*(2*(2*(2*(2*(i1-1)+i2-1)+ # i3-1)+i4-1)+i5-1)+i6-1)+i7-1)+i8-1)+i9 if(ip.eq.i) then if(i1.eq.1) then bl(i,1)= 0.d0 bu(i,1)= prx0n(1) else bl(i,1)= prx0n(1) bu(i,1)= 1.d0 endif if(i2.eq.1) then bl(i,2)= 0.d0 bu(i,2)= prx0n(2) else bl(i,2)= prx0n(2) bu(i,2)= 1.d0 endif if(i3.eq.1) then bl(i,3)= 0.d0 bu(i,3)= prx0n(3) else bl(i,3)= prx0n(3) bu(i,3)= 1.d0 endif if(i4.eq.1) then bl(i,4)= 0.d0 bu(i,4)= prx0n(4) else bl(i,4)= prx0n(4) bu(i,4)= 1.d0 endif if(i5.eq.1) then bl(i,5)= 0.d0 bu(i,5)= prx0n(5) else bl(i,5)= prx0n(5) bu(i,5)= 1.d0 endif if(i6.eq.1) then bl(i,6)= 0.d0 bu(i,6)= prx0n(6) else bl(i,6)= prx0n(6) bu(i,6)= 1.d0 endif if(i7.eq.1) then bl(i,7)= 0.d0 bu(i,7)= prx0n(7) else bl(i,7)= prx0n(7) bu(i,7)= 1.d0 endif if(i8.eq.1) then bl(i,8)= 0.d0 bu(i,8)= prx0n(8) else bl(i,8)= prx0n(8) bu(i,8)= 1.d0 endif if(i9.eq.1) then bl(i,9)= 0.d0 bu(i,9)= prx0n(9) else bl(i,9)= prx0n(9) bu(i,9)= 1.d0 endif endif enddo enddo enddo enddo enddo enddo enddo enddo enddo enddo * do ij=1,nint do i=1,npos vol(i,ij)= 1.d0 do j=1,9 vol(i,ij)= vol(i,ij)*(bu(i,j)-bl(i,j)) enddo vol(i,ij)= vol(i,ij)*xsmx(i,ij) enddo enddo do ij=1,nint volt(ij)= 0.d0 do i=1,npos volt(ij)= volt(ij)+vol(i,ij) enddo enddo * m1= 1 m2= npos ord= 'd' do ij=1,nint do i=1,npos avol(i)= vol(i,ij) enddo ifail= 0 call m01daf(avol,m1,m2,ord,ir,ifail) do i=1,npos vol(ir(i),ij)= avol(i) oxsmx(ir(i),ij)= xsmx(i,ij) do j=1,9 gbl(ir(i),j,ij)= bl(i,j) gbu(ir(i),j,ij)= bu(i,j) enddo enddo enddo * return end * *------------------------------------------------------------------ * function ipack(shat,sm,sp,su,sd,sf) implicit real*8 (a-h,o-z) * common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi * zl= (zm-25.d0)*(zm-25.d0)/shat zu= (zm+25.d0)*(zm+25.d0)/shat bc= 2.5d3/shat sg= 1.d0-sm-sp-su-sd-sf * if(zu.gt.sm.and.sm.gt.zl.and.sp.gt.bc) then ipack= 1 else if(zu.gt.sp.and.sp.gt.zl.and.sm.gt.bc) then ipack= 1 else if(zu.gt.sf.and.sf.gt.zl.and.sg.gt.bc) then ipack= 1 else if(zu.gt.sg.and.sg.gt.zl.and.sf.gt.bc) then ipack= 1 else if(zu.gt.su.and.su.gt.zl.and.sd.gt.bc) then ipack= 1 else if(zu.gt.sd.and.sd.gt.zl.and.su.gt.bc) then ipack= 1 else ipack= 0 endif * return end * *------------------------------------------------------------------ * function ievs(a,b,s1,s2,s3,s4,s5,s6) implicit real*8 (a-h,o-z) * common/wthiggs/hm * c= hm-10.d0 d= hm+10.d0 * if((b.gt.s1.and.s1.gt.a).and.(s2.lt.c.or.s2.gt.d).and. # (s3.lt.c.or.s3.gt.d).and.(s4.lt.c.or.s4.gt.d).and. # (s5.lt.c.or.s5.gt.d).and.(s6.lt.c.or.s6.gt.d)) then ievs= 1 else if((b.gt.s2.and.s2.gt.a).and.(s1.lt.c.or.s1.gt.d).and. # (s3.lt.c.or.s3.gt.d).and.(s4.lt.c.or.s4.gt.d).and. # (s5.lt.c.or.s5.gt.d).and.(s6.lt.c.or.s6.gt.d)) then ievs= 1 else if((b.gt.s3.and.s3.gt.a).and.(s1.lt.c.or.s1.gt.d).and. # (s2.lt.c.or.s2.gt.d).and.(s4.lt.c.or.s4.gt.d).and. # (s5.lt.c.or.s5.gt.d).and.(s6.lt.c.or.s6.gt.d)) then ievs= 1 else if((b.gt.s4.and.s4.gt.a).and.(s1.lt.c.or.s1.gt.d).and. # (s2.lt.c.or.s2.gt.d).and.(s3.lt.c.or.s3.gt.d).and. # (s5.lt.c.or.s5.gt.d).and.(s6.lt.c.or.s6.gt.d)) then ievs= 1 else if((b.gt.s5.and.s5.gt.a).and.(s1.lt.c.or.s1.gt.d).and. # (s2.lt.c.or.s2.gt.d).and.(s3.lt.c.or.s3.gt.d).and. # (s4.lt.c.or.s4.gt.d).and.(s6.lt.c.or.s6.gt.d)) then ievs= 1 else if((b.gt.s6.and.s6.gt.a).and.(s1.lt.c.or.s1.gt.d).and. # (s2.lt.c.or.s2.gt.d).and.(s3.lt.c.or.s3.gt.d).and. # (s4.lt.c.or.s4.gt.d).and.(s5.lt.c.or.s5.gt.d)) then ievs= 1 else ievs= 0 endif * return end * *------------------------------------------------------------------- * subroutine wtopole(p2r,p2i,p3q,fz,fw) implicit real*8(a-h,o-z) character*1 rio * common/wtai/rio common/wtparam/eps,ddelta common/wtfmass2/em2,amm2,tm2,rnm2,uqm2,dqm2,cqm2,sqm2,bqm2, # tqm2,dmy2 * dimension bfl(2),bfh(2) dimension p3q(2),fz(2),fw(2) dimension b0l(2),b1l(2),b21l(2) dimension b0h(2),b1h(2),b21h(2) dimension b0c(2),b1c(2),b21c(2) * call wtosbff(p2r,p2i,dmy2,dmy2,b0l,b1l,b21l) call wtosbff(p2r,p2i,tqm2,tqm2,b0h,b1h,b21h) call wtosbff(p2r,p2i,tqm2,dmy2,b0c,b1c,b21c) * do i=1,2 bfl(i)= 2.d0*b21l(i)-b0l(i) bfh(i)= 2.d0*b21h(i)-b0h(i) enddo bfsr= bfh(1)-bfl(1) bfsi= bfh(2)-bfl(2) * if(rio.eq.'a') then do i=1,2 p3q(i)= 10.d0*bfl(i)+2.d0*bfh(i) enddo else if(rio.eq.'i') then p3q(1)= 0.d0 p3q(2)= 10.d0*bfl(2)+2.d0*bfh(2) endif * if(rio.eq.'a') then fz(1)= -0.5d0*(p2r*bfsr-p2i*bfsi)-3.d0/2.d0*tqm2*b0h(1) else if(rio.eq.'i') then fz(1)= 0.d0 endif fz(2)= -0.5d0*(p2r*bfsi+p2i*bfsr)-3.d0/2.d0*tqm2*b0h(2) * auxr= 3.d0*(b21c(1)+b1c(1))-2.d0*b21h(1)+b0h(1)-b21l(1)- # 3.d0*b1l(1)-b0l(1) auxi= 3.d0*(b21c(2)+b1c(2))-2.d0*b21h(2)+b0h(2)-b21l(2)- # 3.d0*b1l(2)-b0l(2) * if(rio.eq.'a') then fw(1)= 2.d0*(p2r*auxr-p2i*auxi)-3.d0*tqm2*(b1c(1)+b0c(1)) else if(rio.eq.'i') then fw(1)= 0.d0 endif fw(2)= 2.d0*(p2r*auxi+p2i*auxr)-3.d0*tqm2*(b1c(2)+b0c(2)) * return end * *------------------------------------------------------------------------ * subroutine wtobff(p2,rm12,rm22,b0,b1,b21) implicit real*16(a-h,o-z) * common/wtqparam/qpi,qpis,qeps,qdelta * dimension cmp2(2) dimension clnmp2(2) dimension b0(2),b1(2),b21(2) dimension gfpr(3),gfmr(3),gfpi(3),gfmi(3) * cmp2(1)= -p2 cmp2(2)= -qeps * call wtoqroots(p2,rm12,rm22,rpr,rpi,rmr,rmi,omrpr,omrmr) * call wtocqlnx(cmp2,clnmp2) n= 3 call wtocg(n,rpr,rpi,omrpr,gfpr,gfpi) call wtocg(n,rmr,rmi,omrmr,gfmr,gfmi) * auxdel= qdelta-clnmp2(1) b0(1)= auxdel-gfpr(1)-gfmr(1) b0(2)= -clnmp2(2)-gfpi(1)-gfmi(1) b1(1)= -auxdel/2.q0+gfpr(2)+gfmr(2) b1(2)= +clnmp2(2)/2.q0+gfpi(2)+gfmi(2) b21(1)= auxdel/3.q0-gfpr(3)-gfmr(3) b21(2)= -clnmp2(2)/3.q0-gfpi(3)-gfmi(3) * return end * *------------------------------------------------------------------------- * subroutine wtosbff(p2r,p2i,rm12,rm22,b0,b1,b21) implicit real*8(a-h,o-z) * common/wtparam/eps,ddelta common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi * dimension cmp2(2) dimension clnmp2(2) dimension b0(2),b1(2),b21(2) * p2m= p2r*p2r+p2i*p2i * if(rm12.eq.0.d0.and.rm22.eq.0.d0) then if(p2i.eq.0) then cmp2(1)= p2r cmp2(2)= -eps call wtoclnx(cmp2,clnmp2) b0(1)= ddelta+2.d0-clnmp2(1) b0(2)= -clnmp2(2) else O O O O O O O O O O O O O O O O cmp2(1)= 1.d0+p2r/rm12 O cmp2(2)= -eps/rm12% O call wtoclnx(cmp2,clnmp2) O B O b0(1)= ddelta+2.d0-log(rm12)-(1.d0+rm12/p2r)*clnmp2(1)- O b0(2)= -(1.d0+rm12/p2r)*clnmp2(2) else " O cmp2(1)= 1.d0+p2r/rm12 O cmp2(2)= p2i/rm121% O call wtoclnx(cmp2,clnmp2) O  O xr= p2r/p2m( O xi= -p2i/p2m O fctr= 1.d0+rm12*xr O fcti= rm12*xiuF O b0(1)= ddelta+2.d0-log(rm12)-fctr*clnmp2(1)+fcti*clnmp2(2)1 O b0(2)= -fctr*clnmp2(2)-fcti*clnmp2(1)i# O if(p2r.lt.(-rm12)) thend0 O b0(1)= b0(1)+2.d0*pi*rm12*p2i/p2m7 O b0(2)= b0(2)+2.d0*pi*(1.d0+rm12*p2r/p2m)  O endifp O endif1 O else if(rm12.ne.0.d0.and.rm12.eq.rm22) thenl* O a0= rm12*(-ddelta-1.d0+log(rm12)) O if(p2i.eq.0) then O p2d= p2r*p2r+eps*eps' O b2r= 1.d0+4.d0*rm12*p2r/p2d " O b2i= 4.d0*rm12*eps/p2d& O call a02aaf(b2r,b2i,br,bi)C O cmp2(1)= (br*br+bi*bi-1.d0)/((br-1.d0)*(br-1.d0)+bi*bi)n9 O cmp2(2)= -2.d0*bi/((br-1.d0)*(br-1.d0)+bi*bi) % O call wtoclnx(cmp2,clnmp2) 7 O b0(1)= ddelta+2.d0-log(rm12)-(br*clnmp2(1)- O # bi*clnmp2(2))/ O b0(2)= -(br*clnmp2(2)+bi*clnmp2(1))o else  O xr= p2r/p2md O xi= -p2i/p2m" O b2r= 1.d0+4.d0*rm12*xr O b2i= 4.d0*rm12*xi & O call a02aaf(b2r,b2i,br,bi)C O cmp2(1)= (br*br+bi*bi-1.d0)/((br-1.d0)*(br-1.d0)+bi*bi) 9 O cmp2(2)= -2.d0*bi/((br-1.d0)*(br-1.d0)+bi*bi) % O call wtoclnx(cmp2,clnmp2)v7 O b0(1)= ddelta+2.d0-log(rm12)-(br*clnmp2(1)-l O # bi*clnmp2(2))/ O b0(2)= -(br*clnmp2(2)+bi*clnmp2(1)) ( O if(p2r.lt.(-4.d0*rm12)) then& O b0(1)= b0(1)-2.d0*pi*bi& O b0(2)= b0(2)+2.d0*pi*br O endif  O endif O else+ O a01= rm12*(-ddelta-1.d0+log(rm12))p+ O a02= rm22*(-ddelta-1.d0+log(rm22))  O if(p2i.eq.0) then1 O b2r= (p2r-rm22+rm12)*(p2r-rm22+rm12)- ' O # 4.d0*rm12*rm22-eps*eps * O b2i= -2.d0*(p2r+rm22+rm12)*eps& O call a02aaf(b2r,b2i,br,bi)> O cmp2(1)= 0.5d0*(-p2r-rm22-rm12+br)/sqrt(rm22*rm12)3 O cmp2(2)= 0.5d0*(eps+bi)/sqrt(rm22*rm12) % O call wtoclnx(cmp2,clnmp2)-0 O crr= (br*clnmp2(1)-bi*clnmp2(2))/p2r0 O cri= (br*clnmp2(2)+bi*clnmp2(1))/p2rE O b0(1)= ddelta-crr-0.5d0*log(rm22*rm12)+0.5d0*(rm22-rm12)/m* O # p2r*log(rm22/rm12)+2.d0 O b0(2)= -cri O elset O xr= p2r/p2mt O xi= -p2i/p2m O x2r= xr*xr-xi*xi O x2i= 2.d0*xr*xin1 O b2r= (p2r-rm22+rm12)*(p2r-rm22+rm12)-.' O # 4.d0*rm12*rm22-p2i*p2i ) O b2i= 2.d0*(p2r+rm22+rm12)*p2i.& O call a02aaf(b2r,b2i,br,bi)> O cmp2(1)= 0.5d0*(-p2r-rm22-rm12+br)/sqrt(rm22*rm12)4 O cmp2(2)= 0.5d0*(-p2i+bi)/sqrt(rm22*rm12)% O call wtoclnx(cmp2,clnmp2)e0 O crr= (br*clnmp2(1)-bi*clnmp2(2))*xr-/ O # (br*clnmp2(2)+bi*clnmp2(1))*xi-0 O cri= (br*clnmp2(1)-bi*clnmp2(2))*xi+/ O # (br*clnmp2(2)+bi*clnmp2(1))*xrlE O b0(1)= ddelta+crr-0.5d0*log(rm22*rm12)+0.5d0*(rm22-rm12)/0* O # p2r*log(rm22/rm12)+2.d0 O b0(2)= criB O tsti= -(sqrt(rm12)+sqrt(rm22))*(sqrt(rm12)+sqrt(rm22)) O if(p2r.lt.tsti) then1 O b0(1)= b0(1)-2.d0*pi*(br*xi+bi*xr) 1 O b0(2)= b0(2)+2.d0*pi*(br*xr-bi*xi). O endif) O endif O endifs O *t O return O end. O *aF O *--------------------------------------------------------------------- O * = O subroutine wtocff(jflag,op12,op22,os,orm12,orm22,orm32,1% O # oc0,oc1,oc2,oc3)t O implicit real*16(a-h,p-z)4 O implicit real*8(o) O character*1 rio6 O *. O common/wtai/rio=* O common/wtqparam/qpi,qpis,qeps,qdelta O *. O dimension co(2) O ' O dimension tmp3_56(2,2),r3_56(2,2)c+ O dimension b0_12(2),b1_12(2),b21_12(2)t+ O dimension b0_13(2),b1_13(2),b21_13(2)i+ O dimension b0_23(2),b1_23(2),b21_23(2)n- O dimension b22_12(2),b22_13(2),b22_23(2).1 O dimension oc0(2),oc1(2,2),oc2(2,4),oc3(2,6) 0 O dimension r3_34(2,2),r3_13(2,2),r3_42(2,2)6 O dimension tmp3_34(2,2),tmp3_13(2,2),tmp3_42(2,2)8 O dimension r2_32(2,2),tmp2_32(2,2),tmp0(2),tmp24(2)B O dimension xmi(2,2),r1(2,2),tmp1(2,2),r2_13(2,2),tmp2_13(2,2) O *. O qpis= qpi*qpi  O p12= op12*1.d15*1.q-15 O p22= op22*1.d15*1.q-15 O s= os*1.d15*1.q-15 O rm12= orm12*1.d15*1.q-15 O rm22= orm22*1.d15*1.q-15 O rm32= orm32*1.d15*1.q-15 O co(1)= 1.q0p O co(2)= 0.q0 O  O if(jflag.eq.0) then- O rsm= sqrt(p12/s) O  O rsp= sqrt(p22/s) O  O testc= rsm+rspd O if(testc.gt.1.q0) then, O jflag= -1, O return else  O jflag= 1 O endif O endif  O *  O p1dp2= 0.5q0*(s-p12-p22)2 O den= -s*s/4.q0*(1.q0+((p12-p22)*(p12-p22)/s- O # 2.q0*(p12+p22))/s)  O dm12= rm12-rm22, O dm13= rm12-rm32l O dm23= rm22-rm32m O ra= -p22 O rb= -p12 O rc= -2.q0*p1dp2b O rd= dm23+p22 O re= dm12+p12+2.q0*p1dp2b O raa= -p22+dm23 O rbb= p12+dm12  O rcc= -2.q0*p1dp2-p22+dm23  O rdd= -p12+dm12! O discp= +2.q0*sqrt(abs(den))f O discm= -discp  O if(den.lt.0.q0) then5 O arga= (p22-2.q0*(p12+s))*p22+(s-p12)*(s-p12)  O rta= sqrt(arga) O s1a= abs(s/p12) O s2a= abs(s/p22)0 O if(s1a.lt.1.q-10.or.s2a.lt.1.q-10) then O bal= p12-p22+s O if(bal.gt.0.q0) then$ O qal= -0.5q0*(bal+rta) O bep= -s/qal O bem= -qal/p12 O else )$ O qal= -0.5q0*(bal-rta) O bep= -qal/p12 O bem= -s/qal O endif2 O ap= 1.q0-bem O am= 1.q0-bep O qomap= bem O qomam= bep elseb O bal= p12+p22-s O if(bal.gt.0.q0) then$ O qal= -0.5q0*(bal+rta) O ap= -p22/qald O am= -qal/p12  O else $ O qal= -0.5q0*(bal-rta) O ap= -qal/p12+ O am= -p22/qal0 O endife O qomap= 1.q0-ap O qomam= 1.q0-am O endif O if(jflag.eq.1) then O alp= apf O qomalp= qomap O  O disc= discph! O else if(jflag.eq.2) thenp O alp= am  O qomalp= qomam O  O disc= discm elseo O print*,'jflag > 2' O stop O endif O alpi= 0.q01 O y0r= -(rd+re*alp)/disc O  O y0i= 0.q0! O y1r= -(raa+rbb*alp)/disc, O y1i= y0ia# O qomy1r= (rcc+rdd*alp)/disc O  O y2r= y0r/qomalp O y2i= y0i/qomalp O qomy2r= qomy1r/qomalp O y3r= -y0r/alp O y3i= -y0i/alp O qomy3r= y1r/alp O else if(den.gt.0.q0) thenf" O alp= (p12+p22-s)/2.q0/p12% O qomalp= (p12-p22+s)/2.q0/p12l% O calpi= sqrt(den)/p12 ( O if(jflag.eq.1) then O alpi= -calpi O disc= discp O ! O else if(jflag.eq.2) then- O alpi= calpi O disc= discm- else- O print*,'jflag > 2'f O stopb O endif O y0r= -re*alpi/discd O y0i= (rd+re*alp)/disc O y1r= -rbb*alpi/disc) O y1i= (raa+rbb*alp)/disc O  O qomy1r= rdd*alpi/disc* O qomalpm2= qomalp*qomalp+alpi*alpi, O y2r= (y0r*qomalp-y0i*alpi)/qomalpm2, O y2i= (y0r*alpi+y0i*qomalp)/qomalpm22 O qomy2r= (qomy1r*qomalp+y1i*alpi)/qomalpm2! O alpm2= alp*alp+alpi*alpi ' O y3r= -(y0r*alp+y0i*alpi)/alpm2 ( O y3i= -(-y0r*alpi+y0i*alp)/alpm2) O qomy3r= (y1r*alp+y1i*alpi)/alpm2 O  O endif  O *r7 O call wtoqr ts(p12,rm22,rm12,rp1r,rp1i,rm1r,rm1i,,% O # qomrp1r,qomrm1r)05 O call wtoqroots(s,rm32,rm12,rp2r,rp2i,rm2r,rm2i, % O # qomrp2r,qomrm2r)b7 O call wtoqroots(p22,rm32,rm22,rp3r,rp3i,rm3r,rm3i,/% O # qomrp3r,qomrm3r)= O *d O a1= -p12 O a2= -s O a3= -p22 O eb1r= -1.q0-dm12/p12 O eb1i= 0.q0 O ec1i= qeps/p12 O eb2r= -1.q0-dm13/s O eb2i= 0.q0 O ec2i= qeps/s O eb3r= -1.q0-dm23/p22 O eb3i= 0.q0 O ec3i= qeps/p22< O call wtos2(y1r,y1i,qomy1r,a1,eb1r,eb1i,ec1i,rp1r,rp1i,5 O # rm1r,rm1i,qomrp1r,qomrm1r,s21r,s21i) < O call wtos2(y2r,y2i,qomy2r,a2,eb2r,eb2i,ec2i,rp2r,rp2i,5 O # rm2r,rm2i,qomrp2r,qomrm2r,s22r,s22i) < O call wtos2(y3r,y3i,qomy3r,a3,eb3r,eb3i,ec3i,rp3r,rp3i,5 O # rm3r,rm3i,qomrp3r,qomrm3r,s23r,s23i)# O * O *-----scalar 3-point function c0 O *( O if(den.lt.0.q0) then' O tmp0(1)= (s21r-s22r+s23r)/disc ' O tmp0(2)= (s21i-s22i+s23i)/disc  O else if(den.gt.0.q0) thena( O tmp0(1)= (s21i-s22i+s23i)/disc( O tmp0(2)= -(s21r-s22r+s23r)/disc O endif. O *r O f1= dm12-p12 O f2= dm23-s+p12 O * 3 O call wtobff(p12,rm12,rm22,b0_12,b1_12,b21_12)n1 O call wtobff(s,rm12,rm32,b0_13,b1_13,b21_13)m3 O call wtobff(p22,rm22,rm32,b0_23,b1_23,b21_23)  O *  O do i=1,2: O b22_12(i)= -p12*b21_12(i)+0.5q0*(rm12-rm22-p12)*= O # b1_12(i)+0.5q0*co(i)*rm22*(-qdelta-1.q0+x O # log(rm22))f6 O b22_13(i)= -s*b21_13(i)+0.5q0*(rm12-rm32-s)*= O # b1_13(i)+0.5q0*co(i)*rm32*(-qdelta-1.q0+5 O # log(rm32))x: O b22_23(i)= -p22*b21_23(i)+0.5q0*(rm22-rm32-p22)*= O # b1_23(i)+0.5q0*co(i)*rm32*(-qdelta-1.q0+1 O # log(rm32))r O enddo  O *  O xmi(1,1)= p22/den) O xmi(1,2)= -p1dp2/den O xmi(2,1)= xmi(1,2) O xmi(2,2)= p12/den O  O *  O *-----c1 form factors2 O *i O do i= 1,207 O r1(i,1)= 0.5q0*(f1*tmp0(i)+b0_13(i)-b0_23(i))*7 O r1(i,2)= 0.5q0*(f2*tmp0(i)+b0_12(i)-b0_13(i))2 O enddor O do j= 1,2e O do i= 1,29 O tmp1(j,i)= xmi(i,1)*r1(j,1)+xmi(i,2)*r1(j,2)  O enddo O enddo  O *  O *-----c2 form factors1 O *2* O tmp24(1)= 0.25q0-0.5q0*rm12*tmp0(1)+; O # 0.25q0*(b0_23(1)-f1*tmp1(1,1)-f2*tmp1(1,2))($ O tmp24(2)= -0.5q0*rm12*tmp0(2)+; O # 0.25q0*(b0_23(2)-f1*tmp1(2,1)-f2*tmp1(2,2))b O do i= 1,2 E O r2_13(i,1)= 0.5q0*(f1*tmp1(i,1)+b1_13(i)+b0_23(i))-tmp24(i)(< O r2_13(i,2)= 0.5q0*(f2*tmp1(i,1)+b1_12(i)-b1_13(i)) O enddol O do j= 1,2( O do i= 1,2B O tmp2_13(j,i)= xmi(i,1)*r2_13(j,1)+xmi(i,2)*r2_13(j,2) O enddo O enddo2 O do i= 1,2 < O r2_32(i,1)= 0.5q0*(f1*tmp1(i,2)+b1_13(i)-b1_23(i))< O r2_32(i,2)= 0.5q0*(f2*tmp1(i,2)-b1_13(i))-tmp24(i) O enddo  O do j= 1,2  O do i= 1,2B O tmp2_32(j,i)= xmi(i,1)*r2_32(j,1)+xmi(i,2)*r2_32(j,2) O enddo O enddo, O *b O *-----c3 form factors= O *r O do i=1,2= O r3_56(i,1)= 0.5q0*(f1*tmp24(i)+b22_13(i)-b22_23(i))r= O r3_56(i,2)= 0.5q0*(f2*tmp24(i)+b22_12(i)-b22_13(i)) O enddo  O do j= 1,20 O do i= 1,2B O tmp3_56(j,i)= xmi(i,1)*r3_56(j,1)+xmi(i,2)*r3_56(j,2) O enddo O enddo  O *  O do i=1,2A O r3_34(i,1)= 0.5q0*(f1*tmp2_13(i,2)+b21_13(i)+b1_23(i))- " O # tmp3_56(i,2)D O r3_34(i,2)= 0.5q0*(f2*tmp2_13(i,2)-b21_13(i))-tmp3_56(i,1) O enddo  O do j= 1,21 O do i= 1,2B O tmp3_34(j,i)= xmi(i,1)*r3_34(j,1)+xmi(i,2)*r3_34(j,2) O enddo O enddo. O *3 O do i=1,2A O r3_13(i,1)= 0.5q0*(f1*tmp2_13(i,1)+b21_13(i)-b0_23(i))- ' O # 2.q0*tmp3_56(i,1) A O r3_13(i,2)= 0.5q0*(f2*tmp2_13(i,1)+b21_12(i)-b21_13(i))f O enddo  O do j= 1,2e O do i= 1,2B O tmp3_13(j,i)= xmi(i,1)*r3_13(j,1)+xmi(i,2)*r3_13(j,2) O enddo O enddo+ O *2 O do i=1,2A O r3_42(i,1)= 0.5q0*(f1*tmp2_32(i,2)+b21_13(i)-b21_23(i))== O r3_42(i,2)= 0.5q0*(f2*tmp2_32(i,2)-b21_13(i))-2.q0*b" O # tmp3_56(i,2) O enddo/ O do j= 1,2 O  O do i= 1,2B O tmp3_42(j,i)= xmi(i,1)*r3_42(j,1)+xmi(i,2)*r3_42(j,2) O enddo O enddol O *2 O if(rio.eq.'a') then  O do i= 1,2 O oc0(i)= tmp0(i) O oc1(i,1)= tmp1(i,1) O oc1(i,2)= tmp1(i,2)! O oc2(i,1)= tmp2_13(i,1)g! O oc2(i,2)= tmp2_32(i,2)=! O oc2(i,3)= tmp2_13(i,2)r O oc2(i,4)= tmp24(i))! O oc3(i,5)= tmp3_56(i,1)2! O oc3(i,6)= tmp3_56(i,2))! O oc3(i,3)= tmp3_34(i,1)d! O oc3(i,4)= tmp3_34(i,2)r! O oc3(i,1)= tmp3_13(i,1) ! O oc3(i,2)= tmp3_42(i,2)d O enddo O else if(rio.eq.'i') then O oc0(1)= 0.d0r O oc1(1,1)= 0.d0  O oc1(1,2)= 0.d0  O oc2(1,1)= 0.d00 O oc2(1,2)= 0.d0 O  O oc2(1,3)= 0.d00 O oc2(1,4)= 0.d0  O oc3(1,5)= 0.d0  O oc3(1,6)= 0.d0  O oc3(1,3)= 0.d0  O oc3(1,4)= 0.d0  O oc3(1,1)= 0.d0( O oc3(1,2)= 0.d0  O oc0(2)= tmp0(2) O oc1(2,1)= tmp1(2,1) O oc1(2,2)= tmp1(2,2)! O oc2(2,1)= tmp2_13(2,1) ! O oc2(2,2)= tmp2_32(2,2)q! O oc2(2,3)= tmp2_13(2,2)= O oc2(2,4)= tmp24(2) ! O oc3(2,5)= tmp3_56(2,1)m! O oc3(2,6)= tmp3_56(2,2))! O oc3(2,3)= tmp3_34(2,1) ! O oc3(2,4)= tmp3_34(2,2) ! O oc3(2,1)= tmp3_13(2,1)m! O oc3(2,2)= tmp3_42(2,2)n endif 2 O *) O return O endd O *tD O *------------------------------------------------------------------- O *(C O subroutine wtoroots(p2,rm12,rm22,rpr,rpi,rmr,rmi,omrpr,omrmr)2 O implicit real*8(a-h,o-z) O *1 O common/wtparam/eps,ddeltam? O common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi1 O * * O if (rm12.eq.rm22) then" O disc2= (p2+4.d0*rm12)*p2 O else? O disc2= (rm22+2.d0*(p2-rm12))*rm22+(p2+rm12)*(p2+rm12)- O endif  O if (disc2.gt.0.d0) then- O disc= sqrt(disc2)2 O rpi= +eps/disc O rmi= -eps/disc O r1= abs(rm22/p2) O r2= abs(rm22/rm12) O * % O * one of the roots very near to 1o O * 0 O if (r1.lt.1.d-10.or.r2.lt.1.d-10) then O a= -p2 O b= p2+rm12-rm222 O c= rm22#! O if (b.gt.0.d0) then-$ O q= -0.5d0*(b+disc) O rmr= 1.d0-c/qd O rpr= 1.d0-q/a  O omrpr= q/a O omrmr= c/q O else 1$ O q= -0.5d0*(b-disc) O rmr= 1.d0-q/ax O rpr= 1.d0-c/q0 O omrpr= c/q O omrmr= q/a O endif3 O else O a= -p2 O b= p2+rm22-rm12e O c= rm12-! O if (b.gt.0.d0) then-$ O q= -0.5d0*(b+disc) O rpr= c/q O rmr= q/a! O omrpr= 1.d0-rpr ! O omrmr= 1.d0-rmr  O else ,$ O q= -0.5d0*(b-disc) O rpr= q/a O rmr= c/q! O omrpr= 1.d0-rpr ! O omrmr= 1.d0-rmrf O endif2 O endif O else O disc= sqrt(-disc2) O a= -p2 O b= p2+rm22-rm12  O rpr= -b/2.d0/a O rmr= -b/2.d0/a O omrpr= 1.d0-rpr  O omrmr= 1.d0-rmr  O rpi= +disc/2.d0/a2 O rmi= -disc/2.d0/a= O endif O return O endm O *eD O *------------------------------------------------------------------- O *+D O subroutine wtoqroots(p2,rm12,rm22,rpr,rpi,rmr,rmi,omrpr,omrmr) O implicit real*16(a-h,o-z)d O *a* O common/wtqparam/qpi,qpis,qeps,qdelta O *  O if (rm12.eq.rm22) then" O disc2= (p2+4.q0*rm12)*p2 O else? O disc2= (rm22+2.q0*(p2-rm12))*rm22+(p2+rm12)*(p2+rm12). O endif  O if (disc2.gt.0.q0) thend O disc= sqrt(disc2)  O rpi= +qeps/disc1 O rmi= -qeps/disc  O r1= abs(rm22/p2) O r2= abs(rm22/rm12) O * O % O * one of the roots very near to 1i O *p0 O if (r1.lt.1.q-10.or.r2.lt.1.q-10) then O a= -p2 O b= p2+rm12-rm22r O c= rm22i! O if (b.gt.0.q0) thenr$ O q= -0.5q0*(b+disc) O rmr= 1.q0-c/q  O rpr= 1.q0-q/a2 O omrpr= q/a O omrmr= c/q O else n$ O q= -0.5q0*(b-disc) O rmr= 1.q0-q/al O rpr= 1.q0-c/q+ O omrpr= c/q O omrmr= q/a O endif  O else O a= -p2 O b= p2+rm22-rm12 O  O c= rm12(! O if (b.gt.0.q0) then=$ O q= -0.5q0*(b+disc) O rpr= c/q O rmr= q/a! O omrpr= 1.q0-rpro! O omrmr= 1.q0-rmr1 O else 1$ O q= -0.5q0*(b-disc) O rpr= q/a O rmr= c/q! O omrpr= 1.q0-rpr ! O omrmr= 1.q0-rmri O endif= O endifr O else O disc= sqrt(-disc2) O a= -p2 O b= p2+rm22-rm12m O rpr= -b/2.q0/a O rmr= -b/2.q0/a O omrpr= 1.q0-rpr. O omrmr= 1.q0-rmrb O rpi= +disc/2.q0/al O rmi= -disc/2.q0/a O endif+ O return O end( O *lD O *------------------------------------------------------------------- O *cC O subroutine wtos2(y0r,y0i,omy0r,a,ebr,ebi,eci,rpr,rpi,rmr,rmi,b( O # omrpr,omrmr,s2r,s2i) O implicit real*16(a-h,o-z)  O *i* O common/wtqparam/qpi,qpis,qeps,qdelta O *m( O dimension carg(2),comarg(2),cln(2) O *m0 O * rp and rm are the roots of a*x^2+b*x+c = 0? O * eb= b/a and ec= c/a are kept to compute im(rp*rm)= im(ec)#/ O * and im{(y0-rp)(y0-rm)}= im{y0^2+eb*y0+ec}2 O *- O a1= a*eci1. O a2= a*(2.q0*y0r*y0i+ebr*y0i+ebi*y0r+eci) O if (a1.lt.0.q0) then O del= 1.q-90( O else O del= -1.q-90 O endif 2 O if (a2.lt.0.q0) then O delp= 1.q-90 O else O delp= -1.q-902 O endif  O * E O * computes the log which occurs in association with eta-functions- O O  O y0m2= y0r*y0r+y0i*y0i2 O carg(1)= y0r/y0m2  O carg(2)= -y0i/y0m2* O comarg(1)= -(omy0r*y0r-y0i*y0i)/y0m2 O comarg(2)= y0i/y0m2 & O call wtocqlnomx(carg,comarg,cln) O **? O * in calling eta only the sign of the arguments is relevant2 O *1 O a11= -rpi) O a12= -rmi  O a1p= eci O a21= y0i-rpi O a22= y0i-rmi+ O a2p= 2.q0*y0r*y0i+ebr*y0i+ebi*y0r+eci2 O a31= -del  O a32= delp0 O a3p= a*(delp-del)(9 O call wtorlog(y0r,y0i,omy0r,rpr,rpi,omrpr,rfpr,rfpi)r9 O call wtorlog(y0r,y0i,omy0r,rmr,rmi,omrmr,rfmr,rfmi)r O eta1= wtoeta(a11,a12,a1p)- O eta2= wtoeta(a21,a22,a2p)b O eta3= wtoeta(a31,a32,a3p)=1 O s2r= rfpr+rfmr-cln(2)*(eta1-eta2-eta3) (1 O s2i= rfpi+rfmi+cln(1)*(eta1-eta2-eta3) d O return O endr O *)< O *----------------------------------------------------------- O * O A O subroutine wtorlog(z0r,z0i,omz0r,z1r,z1i,omz1r,rfunr,rfuni)( O implicit real*16(a-h,o-z)  O *f* O common/wtqparam/qpi,qpis,qeps,qdelta O *-+ O dimension ca1(2),ca2(2),ca3(2),add(2) 5 O dimension cln1(2),cln2(2),cln3(2),ct(10),sn(10)  O *e O z1m2= z1r*z1r+z1i*z1i  O z1m= sqrt(z1m2) O  O z0m2= z0r*z0r+z0i*z0i- O z0m= sqrt(z0m2)- O zrt= z1m/z0m O *  O * |z1| << |z0| O *f O if (zrt.lt.1.q-10) thenm* O ct(1)= (z1r*z0r+z1i*z0i)/z1m/z0m* O sn(1)= (z0r*z1i-z1r*z0i)/z1m/z0m O do k=2,10 0 O ct(k)= ct(1)*ct(k-1)-sn(1)*sn(k-1)0 O sn(k)= sn(1)*ct(k-1)+ct(1)*sn(k-1) O enddo O  O sumr= ct(10) O sumi= sn(10) O do j=9,1,-1_" O sumr= sumr*zrt+ct(j)" O sumi= sumi*zrt+sn(j) O enddo 2 O sumr= sumr*zrt O sumi= sumi*zrt O a1r= 1.q0+sumr O a1i= sumi,8 O a2r= (z0i*z0i-omz0r*z0r)/z0m2*a1r-z0i/z0m2*a1i8 O a2i= (z0i*z0i-omz0r*z0r)/z0m2*a1i+z0i/z0m2*a1r O oma1r= -sumr O oma2r= 1.q0-a2r, O z01r= z0r-z1r  O z01i= z0i-z1i)" O den= z01r*z01r+z01i*z01i O add(1)= 0.q0 O add(2)= 0.q0 O sign= +1.q0- O else O **1 O * if z0r and/or z1r are roots very near to 1  O * z0r-z1r = (1-z1r)-(1-z0r)m O *  O aomz0r= abs(omz0r) O aomz1r= abs(omz1r) O az1i= abs(z1i)8 O if (aomz0r.lt.1.q-10.or.aomz1r.lt.1.q-10) then O z01r= omz1r-omz0r  O else O z01r= z0r-z1r- O endife O z01i= z0i-z1i O * - O * if z0 and z1r are roots very near to 1 1 O * and |z1i| << 1 , z0i = 0 O *-4 O if (aomz0r.lt.1.q-10.and.aomz1r.lt.1.q-10.6 O # and.az1i.lt.1.q-10.and.z0i.eq.0.q0) then$ O den= z01r*z01r+z1i*z1i O a1r= z01r/z0r  O a1i= -z1i/z0rd O oma1r= z1r/z0r O ca3(1)= -a1r O ca3(2)= -a1i% O call wtocqlnx(ca3,cln3) H O add(1)= -qpis/6.q0-0.5q0*(cln3(1)*cln3(1)-cln3(2)*cln3(2)), O add(2)= -cln3(1)*cln3(2)  O sign= -1.q02 O else & O den= z01r*z01r+z01i*z01i* O a1r= (z0r*z01r+z0i*z01i)/den+ O a1i= (-z0r*z01i+z0i*z01r)/den - O oma1r= -(z01i*z1i+z1r*z01r)/denn O add(1)= 0.q0 O add(2)= 0.q0  O sign= +1.q0q O endif ) O a2r= -(omz0r*z01r-z0i*z01i)/dent( O a2i= (omz0r*z01i+z0i*z01r)/den* O oma2r= (omz1r*z01r-z01i*z1i)/den O endif. O *b? O * in calling eta only the sign of the arguments is relevant  O *a O a1= -z1i O a2= -z01i  O api= z1r*z0i-z0r*z1i O apii= -z0i*omz1r+z1i*omz0r- O call wtospence(a1r,a1i,oma1r,sp1r,sp1i)-- O call wtospence(a2r,a2i,oma2r,sp2r,sp2i)  O ca1(1)= a1r  O ca1(2)= a1i( O ca2(1)= a2r  O ca2(2)= a2i O  O call wtocqlnx(ca1,cln1)  O call wtocqlnx(ca2,cln2)o O etai= wtoeta(a1,a2,api)= O etaii= wtoeta(a1,a2,apii) = O rfunr= sign*sp1r-sp2r-etai*cln1(2)+etaii*cln2(2)+add(1) = O rfuni= sign*sp1i-sp2i+etai*cln1(1)-etaii*cln2(1)+add(2)t O return O end  O * G O *----------------------------------------------------------------------  O *n* O real*16 function wtoeta(z1i,z2i,zpi) O implicit real*16(a-h,o-z)  O *y* O common/wtqparam/qpi,qpis,qeps,qdelta O * 0 O * only the sign of the arguments is relevant O * ; O if (z1i.lt.0.q0.and.z2i.lt.0.q0.and.zpi.gt.0.q0) then O  O wtoeta= 2.q0*qpi O return@ O else if (z1i.gt.0.q0.and.z2i.gt.0.q0.and.zpi.lt.0.q0) then O wtoeta= -2.q0*qpi O return O else O wtoeta= 0.q0 O return O endif1 O end  O *lF O *--------------------------------------------------------------------- O * 2 O subroutine wtospence(xr,xi,omxr,cli2r,cli2i) O implicit real*16(a-h,o-z)i O *c* O common/wtqparam/qpi,qpis,qeps,qdelta O * O O dimension b(0:14),bf(0:14)E O dimension x(2),omx(2),y(2),oy(2),omy(2),z(2),omz(2),t(2),omt(2) < O dimension clnx(2),clnomx(2),clnoy(2),clnz(2),clnomz(2)C O dimension add1(2),add2(2),add3(2),par(2),res(2),ct(15),sn(15) O  O *  O x(1)= xr O x(2)= xi O omx(1)= omxr O omx(2)= -xil O if (xr.lt.0.q0) then O y(1)= omxr O y(2)= -xim O sign1= -1.q0 O call wtocqlnx(x,clnx)l' O call wtocqlnomx(x,omx,clnomx)+@ O add1(1)= qpis/6.q0-clnx(1)*clnomx(1)+clnx(2)*clnomx(2)7 O add1(2)= -clnx(1)*clnomx(2)-clnx(2)*clnomx(1) O O else O y(1)= x(1) O y(2)= x(2) O sign1= 1.q0  O add1(1)= 0.q0  O add1(2)= 0.q0m O endifr O omy(1)= 1.q0-y(1)  O omy(2)= -y(2)r O ym2= y(1)*y(1)+y(2)*y(2) O ym= sqrt(ym2), O if (ym.gt.1.q0) then O z(1)= y(1)/ym2 O z(2)= -y(2)/ym2 O  O sign2= -1.q0 O oy(1)= -y(1) O oy(2)= -y(2)! O call wtocqlnx(oy,clnoy)sA O add2(1)= -qpis/6.q0-0.5q0*((clnoy(1))**2-(clnoy(2))**2) O % O add2(2)= -clnoy(1)*clnoy(2) O O else O z(1)= y(1) O z(2)= y(2) O sign2= 1.q0r O add2(1)= 0.q0s O add2(2)= 0.q0o O endifm O omz(1)= 1.q0-z(1)r O omz(2)= -z(2)  O zr= z(1) O if (zr.gt.0.5q0) then  O t(1)= 1.q0-z(1), O t(2)= -z(2)r O omt(1)= 1.q0-t(1)m O omt(2)= -t(2)2 O sign3= -1.q0 O call wtocqlnx(z,clnz)d' O call wtocqlnomx(z,omz,clnomz)-@ O add3(1)= qpis/6.q0-clnz(1)*clnomz(1)+clnz(2)*clnomz(2)7 O add3(2)= -clnz(1)*clnomz(2)-clnz(2)*clnomz(1)+ O else O t(1)= z(1) O t(2)= z(2) O omt(1)= 1.q0-t(1)2 O omt(2)= -t(2)2 O sign3= 1.q0( O add3(1)= 0.q0, O add3(2)= 0.q0f O endif, O call wtocqlnomx(t,omt,par) O b(0)= 1.q0 O b(1)= -1.q0/2.q0 O b(2)= 1.q0/6.q0  O b(4)= -1.q0/30.q0_ O b(6)= 1.q0/42.q0 O b(8)= -1.q0/30.q0  O b(10)= 5.q0/66.q02 O b(12)= -691.q0/2730.q0 O b(14)= 7.q0/6.q0 O fact= 1.q0 O do n=0,14q O bf(n)= b(n)/fact O fact= fact*(n+2.q0)* O enddo- O parr= par(1) O pari= par(2) O parm2= parr*parr+pari*pari O parm= sqrt(parm2)* O ct(1)= parr/parm O sn(1)= pari/parm O do n=2,15 O , O ct(n)= ct(1)*ct(n-1)-sn(1)*sn(n-1), O sn(n)= sn(1)*ct(n-1)+ct(1)*sn(n-1) O enddo O  O * (@ O res(1)= -((((((((bf(14)*ct(15)*parm2+bf(12)*ct(13))*parm2+? O # bf(10)*ct(11))*parm2+bf(8)*ct(9))*parm2+_= O # bf(6)*ct(7))*parm2+bf(4)*ct(5))*parm2+-A O # bf(2)*ct(3))*(-parm)+bf(1)*ct(2))*(-parm)+ ) O # bf(0)*ct(1))*parm 2@ O res(2)= -((((((((bf(14)*sn(15)*parm2+bf(12)*sn(13))*parm2+? O # bf(10)*sn(11))*parm2+bf(8)*sn(9))*parm2+ = O # bf(6)*sn(7))*parm2+bf(4)*sn(5))*parm2+ A O # bf(2)*sn(3))*(-parm)+bf(1)*sn(2))*(-parm)+-) O # bf(0)*sn(1))*parm A O cli2r= sign1*(sign2*(sign3*res(1)+add3(1))+add2(1))+add1(1) A O cli2i= sign1*(sign2*(sign3*res(2)+add3(2))+add2(2))+add1(2)l O * O return O endd O *=E O *--------------------------------------------------------------------) O * 3 O subroutine wtopoleg(p2r,p2i,pggf,pgglq,pggnp)  O implicit real*8(a-h,o-z) O character*1 rio  O *  O common/wtai/riof O common/wtparam/eps,ddelta ? O common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwipA O common/wtfmass2/em2,amm2,tm2,rnm2,uqm2,dqm2,cqm2,sqm2,bqm2, O  O # tqm2,dmy2  O *t O dimension bfl(2),bfh(2) O dimension pggf(2),pgglq(2)% O dimension b0l(2),b1l(2),b21l(2)(% O dimension b0h(2),b1h(2),b21h(2)1 O *)2 O call wtosbff(p2r,p2i,dmy2,dmy2,b0l,b1l,b21l)2 O call wtosbff(p2r,p2i,tqm2,tqm2,b0h,b1h,b21h) O *i O do i=1,2% O bfl(i)= 2.d0*b21l(i)-b0l(i) % O bfh(i)= 2.d0*b21h(i)-b0h(i)f O enddo) O *1 O if(rio.eq.'a') then  O do i=1,232 O pggf(i)= 12.d0*bfl(i)+16.d0/3.d0*bfh(i)& O pgglq(i)= 44.d0/3.d0*bfl(i) O enddo O else if(rio.eq.'i') then O pggf(1)= 0.d0 O pgglq(1)= 0.d0*0 O pggf(2)= 12.d0*bfl(2)+16.d0/3.d0*bfh(2)$ O pgglq(2)= 44.d0/3.d0*bfl(2) O endif( O *) O if(rio.eq.'i') then  O pggnp= 0.d0 O else if(rio.eq.'a') then O zm2= zm*zmm O p2z= -zm23 O call wtorbff(p2z,em2,em2,b0eez,b1eez,b21eez) O 5 O call wtorbff(p2z,amm2,amm2,b0mmz,b1mmz,b21mmz)*6 O call wtorbff(p2z,tm2,tm2,b0tauz,b1tauz,b21tauz)0 O call wtorbff0(em2,em2,b0ee0,b1ee0,b21ee0)2 O call wtorbff0(amm2,amm2,b0mm0,b1mm0,b21mm0)3 O call wtorbff0(tm2,tm2,b0tau0,b1tau0,b21tau0)b O bfez= 2.d0*b21eez-b0eez O bfmz= 2.d0*b21mmz-b0mmz" O bftauz= 2.d0*b21tauz-b0tauz O bfe0= 2.d0*b21ee0-b0ee0 O bfm0= 2.d0*b21mm0-b0mm0" O bftau0= 2.d0*b21tau0-b0tau0$ O pggl= 4.d0*(bfez+bfmz+bftauz)% O pggl0= 4.d0*(bfe0+bfm0+bftau0)o O b= 0.00299d0  O c= 1.d0? O an= -b*log(1.d0+c*zm2)-1.d0/alphai/4.d0/pi*(pggl-pggl0)+i( O # 1.d0-128.89d0/alphai  O ap2x= abs(p2r)  O px= sqrt(ap2x)  O if(px.lt.0.3d0) then1 O a= 0.d0i O b= 0.00835d0 O c= 1.d0p O else if(px.lt.3.d0) thent O a= 0.d0  O b= 0.00238d0 O c= 3.927d0! O else if(px.lt.100.d0) then  O a= 0.00165d0 O b= 0.00299d0 O c= 1.d0o! O else if(px.gt.100.d0) then= O a= 0.00221d0 O b= 0.00293d0 O c= 1.d0 O endif O *0 O da= an-0.00165q0= O a= a+da3 O pggnp= 4.d0*pi*alphai*(a+b*log(1.d0+c*ap2x)) O endif= O *d O return O end0 O * J O *------------------------------------------------------------------------- O *p% O subroutine wtopolez0(pggf0,fw0)  O implicit real*8(a-h,o-z) O character*1 rio_ O *2 O common/wtai/rio2 O common/wtparam/eps,ddelta2? O common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi3A O common/wtfmass2/em2,amm2,tm2,rnm2,uqm2,dqm2,cqm2,sqm2,bqm2,  O # tqm2,dmy2  O *  O if(rio.eq.'i') then  O pggf0= 0.d0p O fw0= 0.d0  O return O else if(rio.eq.'a') then0 O call wtorbff0(em2,em2,b0ee0,b1ee0,b21ee0)2 O call wtorbff0(amm2,amm2,b0mm0,b1mm0,b21mm0)3 O call wtorbff0(tm2,tm2,b0tau0,b1tau0,b21tau0)r2 O call wtorbff0(tqm2,tqm2,b0tt0,b1tt0,b21tt0) O *  O bfe= 2.d0*b21ee0-b0ee0 O  O bfm= 2.d0*b21mm0-b0mm0,! O bftau= 2.d0*b21tau0-b0tau0  O bft= 2.d0*b21tt0-b0tt0  O *n1 O pggf0= 4.d0*(bfe+bfm+bftau)+16.d0/3.d0*bft  O * 2 O call wtorbff0(tqm2,dmy2,b0tb0,b1tb0,b21tb0) O *2$ O fw0= -3.d0*tqm2*(b0tb0+b1tb0) return O endif  O *t O end  O * > O *------------------------------------------------------------- O *2* O subroutine wtofst(nt,st,fvect,iflag) O implicit real*8(a-h,o-z) O * O ? O common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi A O common/wtfmass2/em2,amm2,tm2,rnm2,uqm2,dqm2,cqm2,sqm2,bqm2,  O # tqm2,dmy2- O *  O dimension zs(2),zl(2)d O dimension st(nt),fvect(nt)@ O dimension pggfz(2),pgglqz(2),p3qz(2),fzz(2),p3qw(2),fw(2), O # fzw(2),fwz(2)q O *  O tqm2= st(1) O  O *  O *-----alpha(zm)- O *d O wm2= wm*wm O zm2= zm*zm O pzr= -zm2  O pzi= 0.d0d0 O call wtopoleg(pzr,pzi,pggfz,pgglqz,pggnpz) O call wtopolez0(pggf0,fw0)n1 O x= alphai-0.25d0*(pggfz(1)-pggf0+pggnpz)/pi ( O y= -0.25d0*(pggfz(2)+pgglqz(2))/pi O alzr= x/(x*x+y*y). O alzi= -y/(x*x+y*y) O *  O *-----re(1/g^2(zm))) O *  O apis= 16.d0*pisq O p2zr= -zm2 O p2zi= 0.d0* O call wtopole(p2zr,p2zi,p3qz,fzz,fwz)0 O bx= 1.d0/8.d0/gf/zm2+(fw0-fzz(1))/zm2/apis O xih= -p3qz(2)/apis O ac= 4.d0*pi*alzr O bc= -1.d0-8.d0*pi*xih*alzi" O cc= -4.d0*pi*xih*xih*alzr+bx O ifail= 0' O call c02ajf(ac,bc,cc,zs,zl,ifail)i O g2z= zs(1) O * O  O p2r= -wm2  O p2i= 0.d0i' O call wtopole(p2r,p2i,p3qw,fzw,fw)m O *r@ O fvect(1)= 1.d0+0.5d0*gf*wm2/pis*((fw0-fw(1))/wm2+p3qw(1))-. O # 8.d0*gf*wm2*(g2z+p3qz(1)/apis) O * O return O end O  O * C O *------------------------------------------------------------------ O  O *-) O subroutine wtofsw(n,sw,fvecw,iflag)- O implicit real*8(a-h,o-z) O *r? O common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi  O *l O dimension sw(n),fvecw(n)8 O dimension p3qw(2),fzw(2),fw(2),p3q(2),fz(2),fcw(2) O *m O wm2= wm*wm O swr= sw(1) O swi= sw(2) O wmu= sqrt(swr) O wlg= -swi/wmu- O *2 O call wtopolez0(pggf0,fw0)  O p2r= -wm2  O p2i= 0.d0)' O call wtopole(p2r,p2i,p3qw,fzw,fw)  O p2r= -swrd O p2i= -swim& O call wtopole(p2r,p2i,p3q,fz,fcw) O * O  O x0= 1.d0/8.d0/gf/wm2@ O xr= 1.d0+0.5d0*gf*wm2/pis*((fw0-fw(1))/wm2+p3qw(1)-p3q(1))" O xi= -0.5d0*gf*wm2/pis*p3q(2) O xr= x0*xr2 O xi= x0*xi= O *+9 O fvecw(1)= gf*(swr*xr-swi*xi+(fcw(1)-fw0)/16.d0/pis)0 O # -1.d0/8.d0. O fvecw(2)= swr*xi+swi*xr+fcw(2)/16.d0/pis O */ O return O end= O *q; O *----------------------------------------------------------= O *q) O subroutine wtofsz(n,sz,fvecz,iflag)= O implicit real*8(a-h,o-z) O *  O common/wtafsz/g2z ( O common/wtiaz/alzr,alzi,alzir,alzii? O common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi  O *  O dimension sz(n),fvecz(n)8 O dimension p3qz(2),fzz(2),p3q(2),fz(2),fwz(2),fw(2): O dimension pggfz(2),pgglqz(2),pggf(2),pgglq(2),evl(2) O *b O szr= sz(1) O szi= sz(2) O *  O zm2= zm*zm O apis= 16.d0*pis  O call wtopolez0(pggf0,fw0)  O pzr= -zm21 O pzi= 0.d0l0 O call wtopoleg(pzr,pzi,pggfz,pgglqz,pggnpz)( O call wtopole(pzr,pzi,p3qz,fzz,fwz) O p2r= -szr  O p2i= -sziq- O call wtopoleg(p2r,p2i,pggf,pgglq,pggnp) % O call wtopole(p2r,p2i,p3q,fz,fw)  O do i=1,25 O evl(i)= pggfz(i)-pggf(i)+pgglqz(i)-pgglq(i)+ O enddo ! O axr= alzir+0.25d0*evl(1)/pim! O axi= alzii+0.25d0*evl(2)/pi- O alr= axr/(axr*axr+axi*axi)! O ali= -axi/(axr*axr+axi*axi)  O * # O xr= g2z+(p3qz(1)-p3q(1))/apis  O xi= -p3q(2)/apis O *-7 O yr= xr-4.d0*pi*(alr*(xr*xr-xi*xi)-2.d0*ali*xr*xi)-7 O yi= xi-4.d0*pi*(2.d0*alr*xr*xi+ali*(xr*xr-xi*xi))e O *e= O fvecz(1)= 8.d0*gf*(szr*yr-szi*yi+(fz(1)-fw0)/apis)-1.d0 O ( O fvecz(2)= szr*yi+szi*yr+fz(2)/apis O *o O return O endp O *dE O *--------------------------------------------------------------------  O * , O subroutine wtocg(n,zr,zi,omzr,gfr,gfi) O implicit real*16(a-h,o-z)r O *m* O common/wtqparam/qpi,qpis,qeps,qdelta O *^ O dimension gfr(n),gfi(n)c3 O dimension ca(2,10),aux(2),zp(2),ct(16),sn(16) B O dimension a(2,4),ra(16),z(2),omz(2),oz(2),clnomz(2),clnoz(2) O * + O data ra/-0.2q+0,0.666666666666667q-1,2> O # -0.952380952380952q-2,-0.396825396825397q-3,> O # 0.317460317460317q-3,-0.132275132275132q-4,= O # -0.962000962000962q-5,0.105218855218855q-5,0> O # 0.266488361726450q-6,-0.488745528428000q-7,= O # -0.675397500794000q-8,0.190720263471000q-8,m? O # 0.153663007690000q-9,-0.679697905790000q-10,e? O # -0.293683556000000q-11,0.228836696000000q-11/p O * O  O z(1)= zr O z(2)= zi O omz(1)= omzr O omz(2)= -zii) O if (zr.eq.0.q0.and.zi.eq.0.q0) thenc O do k=1,n O gfr(k)= -1.q0/k/k O gfi(k)= 0.q0o O enddor O return. O else if (zr.eq.1.q0.and.zi.eq.0.q0) then O a(1,1)= -1.q0 O a(2,1)= qpi O do j=2,4a7 O a(1,j)= ((j-1.q0)*a(1,j-1)-1.q0/j)/jf. O a(2,j)= (j-1.q0)/j*a(2,j-1) O enddot O do k=1,n O gfr(k)= a(1,k) O gfi(k)= a(2,k) O enddo  O return O else z O zmod2= zr*zr+zi*zi O zmod= sqrt(zmod2)f# O call wtocqlnx(omz,clnomz), O *  O * |z| < 4  O * 2 O if (zmod.lt.4.q0) then  O oz(1)= -z(1)1 O oz(2)= -z(2) 1$ O call wtocqlnx(oz,clnoz)F O ca(1,1)= omz(1)*clnomz(1)-omz(2)*clnomz(2)+z(1)*clnoz(1)-( O # z(2)*clnoz(2)-1.q0F O ca(2,1)= omz(1)*clnomz(2)+omz(2)*clnomz(1)+z(1)*clno