diff --git a/Looptools/B/ffcb0.F b/Looptools/B/ffcb0.F --- a/Looptools/B/ffcb0.F +++ b/Looptools/B/ffcb0.F @@ -1,836 +1,837 @@ #include "externals.h" #include "types.h" * $Id: ffcb0.f,v 1.11 1996/07/18 10:49:04 gj Exp $ *###[ ffcb0: subroutine ffcb0(cb0,cp,cma,cmb,ier) ***#[*comment:*********************************************************** * * * calculates the the two-point function (cf 't Hooft and Veltman) * * we include an overall factor 1/(i*pi^2) relative to FormF * * * * Input: cp (complex) k2, in B&D metric * * cma (complex) mass2, re>0, im<0. * * cmb (complex) mass2, re>0, im<0. * * * * Output: cb0 (complex) B0, the two-point function, * * ier (integer) number of digits lost in calculation * * * * Calls: ffcb0p,ffxb0p * * * ***#]*comment:*********************************************************** * #[ declarations: implicit none * * arguments * integer ier ComplexType cb0,cp,cma,cmb * * local variables * integer init,ithres,i,j,nschsa logical lreal ComplexType cmamb,cmap,cmbp,cm,c,cb0p,cqi(3),cqiqj(3,3) RealType absc,xp,xma,xmb,sprec,smax save init * * common blocks * #include "ff.h" * * statement function * absc(c) = abs(Re(c)) + abs(Im(c)) * * data * data init /0/ * * #] declarations: * #[ the real cases: * if ( Im(cma) .eq. 0 .and. Im(cmb) .eq. 0 .and. + Im(cp).eq.0 ) then lreal = .TRUE. elseif ( nschem.le.4 ) then lreal = .TRUE. if ( init.eq.0 ) then init = 1 print *,'ffcb0: nschem <= 4, ignoring complex masses: ', + nschem endif elseif ( nschem.le.6 ) then if ( init.eq.0 ) then init = 1 print *,'ffcb0: nschem = 5,6 complex masses near ', + 'threshold: ',nschem endif cqi(1) = cma cqi(2) = cmb cqi(3) = cp cqiqj(1,2) = cma - cmb cqiqj(2,1) = -cqiqj(1,2) cqiqj(1,3) = cma - cp cqiqj(3,1) = -cqiqj(1,3) cqiqj(2,3) = cmb - cp cqiqj(3,2) = -cqiqj(2,3) cqiqj(1,1) = 0 cqiqj(2,2) = 0 cqiqj(3,3) = 0 call ffthre(ithres,cqi,cqiqj,3,1,2,3) if ( ithres.eq.0 .or. ithres.eq.1 .and. nschem.eq.5 ) then lreal = .TRUE. else lreal = .FALSE. endif else lreal = .FALSE. endif if ( lreal ) then xp = Re(cp) xma = Re(cma) xmb = Re(cmb) sprec = precx precx = precc call ffxb0(cb0,xp,xma,xmb,ier) precx = sprec if ( ldot ) then do 120 j=1,3 do 110 i=1,3 cfpij2(i,j) = fpij2(i,j) 110 continue 120 continue endif return endif * * #] the real cases: * #[ get differences: * cmamb = cma - cmb cmap = cma - cp cmbp = cmb - cp * * #] get differences: * #[ calculations: * * no more schem-checking, please... * nschsa = nschem nschem = 7 call ffcb0p(cb0p,cp,cma,cmb,cmap,cmbp,cmamb,ier) nschem = nschsa if ( cma .eq. 0 ) then if ( cmb .eq. 0 ) then cm = 1 else cm = cmb**2 endif elseif ( cmb .eq. 0 ) then cm = cma**2 else cm = cma*cmb endif if ( mudim .ne. 0 ) cm = cm/Re(mudim)**2 if ( absc(cm) .gt. xclogm ) then cb0 = Re(delta) - cb0p - log(cm)/2 smax = max(abs(delta),absc(cb0p),absc(log(cm))/2) else call fferr(3,ier) cb0 = -cb0p + Re(delta) endif * #] calculations: *###] ffcb0: end *###[ ffcb0p: subroutine ffcb0p(cb0p,cp,cma,cmb,cmap,cmbp,cmamb,ier) ***#[*comment:*********************************************************** * * * calculates the main part of the two-point function (cf 't * * Hooft and Veltman) for all possible cases: masses equal, * * unequal, equal to zero, real or complex (with a negative * * imaginary part). I think it works. * * Has been checked against FormF for all parameter space. * * Only problems with underflow for extreme cases. VERY OLD CODE. * * * * Input: cp (complex) k2, in B&D metric * * cma (complex) mass2, re>0, im<0. * * cmb (complex) mass2, re>0, im<0. * * cmap/b (complex) cma/b - cp * * cmamb (complex) cma - cmb * * * * Output: cb0p (complex) B0, the two-point function, * * minus log(cm/mu), delta and the * * factor -ipi^2. * * ier (integer) 0=ok, 1=numerical problems, 2=error * * * * Calls: (z/a)log, atan. * * * ***#]*comment:*********************************************************** * #[ declarations: implicit none * * arguments * integer ier ComplexType cb0p,cp,cma,cmb,cmap,cmbp,cmamb * * local variables * integer i,j,initeq,initn1,n1,n2,nffeta,nffet1,init, + ithres,is1 logical lreal RealType xp,ax,ay,ffbnd, + xprceq,bdeq01,bdeq05,bdeq11,bdeq17,bdeq25, + xprcn1,bdn101,bdn105,bdn110,bdn115, + xprnn2,bdn201,bdn205,bdn210,bdn215, + xpneq(30),xpnn1(30), + absc,sprec,xma,xmb,dmap,dmbp,dmamb,smax ComplexType cm,cmp,cm1,cm2,cm1m2, + cm1p,cm2p,cs,cs1,cs2,cx,cy,csom,clam,cslam,clogmm, + zfflo1,c,zm,zp,zm1,zp1,zfflog,cqi(3), + cqiqj(3,3),cpiDpj(3,3),ck,clamr,cslamr,zmr,zpr,zm1r,zp1r save initeq,initn1,xpneq,xpnn1,init, + xprceq,bdeq01,bdeq05,bdeq11,bdeq17,bdeq25, + xprcn1,bdn101,bdn105,bdn110,bdn115, + xprnn2,bdn201,bdn205,bdn210,bdn215 *FOR ABSOFT ONLY * ComplexType csqrt * external csqrt * * common blocks * #include "ff.h" * * data * data xprceq /-1./ data xprcn1 /-1./ data xprnn2 /-1./ data initeq /0/ data initn1 /0/ data init /0/ * * statement function * absc(c) = abs(Re(c)) + abs(Im(c)) * * #] declarations: * #[ fill some dotproducts: * call ffcot2(cpiDpj,cp,cma,cmb,cmap,cmbp,cmamb,ier) if ( ldot ) then do 20 i=1,3 do 10 j=1,3 cfpij2(j,i) = cpiDpj(j,i) fpij2(j,i) = Re(cpiDpj(j,i)) 10 continue 20 continue endif * * #] fill some dotproducts: * #[ the real cases: * if ( Im(cma) .eq. 0 .and. Im(cmb) .eq. 0 .and. + Im(cp).eq.0 ) then lreal = .TRUE. elseif ( nschem.le.4 ) then lreal = .TRUE. if( init.eq.0 ) then init = 1 print *,'ffcb0p: nschem <= 4, ignoring complex masses:', + nschem endif elseif ( nschem.le.6 ) then if( init.eq.0 ) then init = 1 print *,'ffcb0p: nschem = 4,6 complex masses near ', + 'threshold: ',nschem endif cqi(1) = cma cqi(2) = cmb cqi(3) = cp cqiqj(1,2) = cmamb cqiqj(2,1) = -cqiqj(1,2) cqiqj(1,3) = cmap cqiqj(3,1) = -cqiqj(1,3) cqiqj(2,3) = cmbp cqiqj(3,2) = -cqiqj(2,3) cqiqj(1,1) = 0 cqiqj(2,2) = 0 cqiqj(3,3) = 0 call ffthre(ithres,cqi,cqiqj,3,1,2,3) if ( ithres.eq.0 .or. ithres.eq.1 .and. nschem.eq.5 ) then lreal = .TRUE. else lreal = .FALSE. endif else lreal = .FALSE. endif if ( lreal ) then xp = Re(cp) xma = Re(cma) xmb = Re(cmb) dmap = Re(cmap) dmbp = Re(cmbp) dmamb = Re(cmamb) sprec = precx precx = precc call ffxb0p(cb0p,xp,xma,xmb,dmap,dmbp,dmamb,ier) precx = sprec if ( ldot ) then do 120 j=1,3 do 110 i=1,3 cfpij2(i,j) = fpij2(i,j) 110 continue 120 continue endif return endif * * #] the real cases: * #[ which case: * * sort according to the type of mass combination encountered: * 200: one equal to zero, 300: both equal, 400: rest. * if ( cma .eq. 0 ) then if ( cmb .eq. 0 ) then goto 100 endif cm = cmb cmp = cmbp goto 200 endif if ( cmb .eq. 0 ) then cm = cma cmp = cmap goto 200 endif if ( cma .eq. cmb ) then cm = cma cmp = cmap goto 300 endif if ( Re(cma) .lt. Re(cmb) ) then cm2 = cma cm1 = cmb cm1m2 = -cmamb cm1p = cmbp cm2p = cmap is1 = 2 else cm1 = cma cm2 = cmb cm1m2 = cmamb cm1p = cmap cm2p = cmbp is1 = 1 endif goto 400 * #] which case: * #[ both masses equal to zero: 100 continue if ( absc(cp) .gt. xclogm ) then if ( Re(cp).gt.0 ) then cb0p = log(cp) - c2ipi/2 - 2 else cb0p = log(-cp) - 2 endif else cb0p = 0 call fferr(7,ier) endif return * #] both masses equal to zero: * #[ one mass zero: 200 continue * * special case cp = 0, checked 25-oct-1991 * if ( cp .eq. 0 ) then cb0p = -1 goto 990 endif * * Normal case: * cs1 = cp/cm cs2 = cmp/cm * make sure we get the right Riemann sheet! if ( absc(cs1) .lt. xloss ) then cs = zfflo1(cs1,ier) elseif ( Re(cs2).gt.0 ) then cs = zfflog(cs2,0,czero,ier) else cs = zfflog(-cs2,0,czero,ier) cs = cs - c2ipi/2 endif cs = -cs*cmp/cp cb0p = cs - 2 goto 990 * #] one mass zero: * #[ both masses equal: 300 continue * * Both masses are equal. Not only this speeds up things, some * cancellations have to be avoided as well. Checked 25-oct-1991. * -#[ taylor expansion: * * first this special case * if ( absc(cp) .lt. 8*xloss*absc(cm) ) then * * a Taylor expansion seems appropriate as the result will go * as k^2 but seems to go as 1/k !! * * #[ data and bounds: if ( initeq .eq. 0 ) then initeq = 1 xpneq(1) = 1/6D0 do 1 i=2,30 xpneq(i) = xpneq(i-1)*Re(i-1)/Re(2*(2*i+1)) 1 continue endif if (xprceq .ne. precc ) then * * calculate the boundaries for the number of terms to be * included in the taylorexpansion * xprceq = precc sprec = precx precx = precc bdeq01 = ffbnd(1,1,xpneq) bdeq05 = ffbnd(1,5,xpneq) bdeq11 = ffbnd(1,11,xpneq) bdeq17 = ffbnd(1,17,xpneq) bdeq25 = ffbnd(1,25,xpneq) precx = sprec endif * #] data and bounds: cx = cp/cm ax = absc(cx) if ( ax .gt. bdeq17 ) then csom = cx*(Re(xpneq(18)) + cx*(Re(xpneq(19)) + + cx*(Re(xpneq(20)) + cx*(Re(xpneq(21)) + + cx*(Re(xpneq(22)) + cx*(Re(xpneq(23)) + + cx*(Re(xpneq(24)) + cx*(Re(xpneq(25)) )))))))) else csom = 0 endif if ( ax .gt. bdeq11 ) then csom = cx*(Re(xpneq(12)) + cx*(Re(xpneq(13)) + + cx*(Re(xpneq(14)) + cx*(Re(xpneq(15)) + + cx*(Re(xpneq(16)) + cx*(Re(xpneq(17)) + csom )))) + )) endif if ( ax .gt. bdeq05 ) then csom = cx*(Re(xpneq(6)) + cx*(Re(xpneq(7)) + + cx*(Re(xpneq(8)) + cx*(Re(xpneq(9)) + + cx*(Re(xpneq(10)) + cx*(Re(xpneq(11)) + csom )))))) endif if ( ax .gt. bdeq01 ) then csom = cx*(Re(xpneq(2)) + cx*(Re(xpneq(3)) + + cx*(Re(xpneq(4)) + cx*(Re(xpneq(5)) + csom )))) endif cb0p = -cx*(Re(xpneq(1))+csom) goto 990 endif * -#] taylor expansion: * -#[ normal case: * * normal case. first determine if the arguments of the logarithm * has positive real part: (we assume Re(cm) > Im(cm) ) * call ffclmb(clam,-cp,-cm,-cm,cmp,cmp,czero) cslam = sqrt(clam) call ffcoot(zm,zp,cone,chalf,cm/cp,cslam/(2*cp),ier) cs1 = zp/zm if ( absc(cs1-1) .lt. xloss ) then * In this case a quicker and more accurate way is to * calculate log(1-cx). cs2 = cp - cslam if ( absc(cs2) .lt. xloss*absc(cp) ) then cs2 = -cslam*(cp+cslam)/(4*cp*cm) else cs2 = -2*cslam/cs2 endif cs = zfflo1(cs2/(2*cm),ier) else * finally the normal case cs = zfflog(cs1,0,czero,ier) endif cs = cslam*cs/cp cb0p = cs - 2 * * eta terms * n1 = nffet1(zp,1/zm,cs1,ier) if ( Im(cp).eq.0 ) then n2 = nffet1(-zp,-1/zm,cs1,ier) else * use the onshell expression to get the correct continuation ck = Re(cp) call ffclmb(clamr,-ck,-cm,-cm,cm-ck,cm-ck,czero) cslamr = sqrt(clamr) call ffcoot(zmr,zpr,cone,chalf,cm/ck,cslamr/(2*ck),ier) if ( absc(zm-zmr)+absc(zp-zpr).gt.absc(zm-zpr)+absc(zp-zmr) + ) then cs1 = zmr zmr = zpr zpr = cs1 endif if ( Im(zmr).eq.0 .or. Im(zpr).eq.0 ) then if ( Re(zpr).gt.Re(zmr) ) then n2 = +1 else n2 = -1 endif else n2 = nffeta(-zpr,-1/zmr,ier) endif endif if ( n1+n2 .ne. 0 ) + cb0p = cb0p - cslam*c2ipi*(n1+n2)/(2*cp) * also superfluous - just to make sure goto 990 * -#] normal case: * * #] both masses equal: * #[ unequal nonzero masses: 400 continue * -#[ get log(xm2/xm1): cx = cm2/cm1 c = cx-1 if ( 1/absc(cx) .lt. xclogm ) then call fferr(6,ier) clogmm = 0 elseif ( absc(c) .lt. xloss ) then clogmm = zfflo1(cm1m2/cm1,ier) else clogmm = log(cx) endif * -#] get log(xm2/xm1): * -#[ cp = 0: * * first a special case * if ( cp .eq. 0 ) then cs2 = ((cm2+cm1) / cm1m2)*clogmm * save the factor 1/2 for the end cs = - cs2 - 2 if ( absc(cs) .lt. xloss*2 ) then * Taylor expansions: choose which one cx = cm1m2/cm1 ax = absc(cx) if ( ax .lt. .15 .or. precc .gt. 1.E-8 .and. ax + .lt. .3 ) then * #[ taylor 1: * * This is the simple Taylor expansion 'n1' * *--#[ data and bounds: * get the coefficients of the taylor expansion if ( initn1 .eq. 0 ) then initn1 = 1 do 410 i = 1,30 - 410 xpnn1(i)=Re(i)/Re((i+1)*(i+2)) + xpnn1(i)=Re(i)/Re((i+1)*(i+2)) + 410 continue endif * determine the boundaries for 1,5,10,15 terms if ( xprcn1 .ne. precc ) then xprcn1 = precc sprec = precx precx = precc bdn101 = ffbnd(1,1,xpnn1) bdn105 = ffbnd(1,5,xpnn1) bdn110 = ffbnd(1,10,xpnn1) bdn115 = ffbnd(1,15,xpnn1) precx = sprec endif *--#] data and bounds: * calculate: if ( ax .gt. bdn110 ) then cs = cx*(Re(xpnn1(11)) + cx*(Re(xpnn1(12)) + + cx*(Re(xpnn1(13)) + cx*(Re(xpnn1(14)) + + cx*(Re(xpnn1(15))) )))) else cs = 0 endif if ( ax .gt. bdn105 ) then cs = cx*(Re(xpnn1(6)) + cx*(Re(xpnn1(7)) + + cx*(Re(xpnn1(8)) + cx*(Re(xpnn1(9)) + + cx*(Re(xpnn1(10)) + cs))))) endif if ( ax .gt. bdn101 ) then cs = cx*(Re(xpnn1(2)) + cx*(Re(xpnn1(3)) + + cx*(Re(xpnn1(4)) + cx*(Re(xpnn1(5)) + + cs)))) endif cs = cx*cx*(Re(xpnn1(1)) + cs) * #] taylor 1: else * #[ taylor 2: * * This is the more complicated exponential Taylor * expansion 'n2' * * #[ bounds: * determine the boundaries for 1,5,10,15 terms for this * Taylor expansion (starting at i=4) * if ( xprnn2 .ne. precc ) then xprnn2 = precc sprec = precx precx = precc bdn201 = ffbnd(4,1,xinfac) bdn205 = ffbnd(4,5,xinfac) bdn210 = ffbnd(4,10,xinfac) bdn215 = ffbnd(4,15,xinfac) precx = sprec endif * #] bounds: * calculate: cy = 2*cx/(2-cx) ay = absc(cy) if ( ay .gt. bdn210 ) then cs = cy*(Re(xinfac(14)) + cy*(Re(xinfac(15)) + + cy*(Re(xinfac(16)) + cy*(Re(xinfac(17)) + + cy*(Re(xinfac(18))))))) else cs = 0 endif if ( ay .gt. bdn205 ) then cs = cy*(Re(xinfac(9)) + cy*(Re(xinfac(10)) + + cy*(Re(xinfac(11)) + cy*(Re(xinfac(12)) + + cy*(Re(xinfac(13)) + cs))))) endif if ( ay .gt. bdn201 ) then cs = cy*(Re(xinfac(5)) + cy*(Re(xinfac(6)) + + cy*(Re(xinfac(7)) + cy*(Re(xinfac(8)) + + cs)))) endif cs = (1-cx)*cy**4 * (Re(xinfac(4)) + cs) cs = cx*cy**2*(1+cy)/12 - cs cs = - 2*zfflo1(cs,ier)/cy * #] taylor 2: endif endif cb0p = cs/2 goto 990 endif * -#] cp = 0: * -#[ normal case: * * (programmed anew 28-oct-1991) * call ffclmb(clam,cm1,cm2,cp,cm1m2,cm1p,cm2p) cslam = sqrt(clam) if ( is1.eq.1 ) then cs = +cpiDpj(2,3) else cs = -cpiDpj(1,3) endif call ffcoot(zm,zp,cp,cs,cm2,cslam/2,ier) zm1 = 1-zm zp1 = 1-zp if ( absc(zm1) .lt. xloss .or. absc(zp1) .lt. xloss ) then if ( is1.eq.1 ) then cs = -cpiDpj(1,3) else cs = +cpiDpj(2,3) endif call ffcoot(zp1,zm1,cp,cs,cm1,cslam/2,ier) if ( abs(Im(zm)) .lt. abs(Im(zm1)) ) then zm = ToComplex(Re(zm),-Im(zm1)) else zm1 = ToComplex(Re(zm1),-Im(zm)) endif if ( abs(Im(zp)) .lt. abs(Im(zp1)) ) then zp = ToComplex(Re(zp),-Im(zp1)) else zp1 = ToComplex(Re(zp1),-Im(zp)) endif endif if ( Im(cp).ne.0 ) then * compute roots for Im(cp).eq.0 for continuation terms. ck = Re(cp) call ffclmb(clamr,cm1,cm2,ck,cm1m2,cm1-ck,cm2-ck) cslamr = sqrt(clamr) if ( absc(cslamr-cslam).gt.absc(cslamr+cslam) ) + cslamr = -cslamr cs = (cm2-cm1+ck)/2 call ffcoot(zmr,zpr,ck,cs,cm2,cslamr/2,ier) zm1r = 1-zmr zp1r = 1-zpr if ( absc(zm1r) .lt. xloss .or. absc(zp1r) .lt. xloss ) then cs = -(cm2-cm1-ck)/2 call ffcoot(zp1r,zm1r,ck,cs,cm1,cslamr/2,ier) if ( abs(Im(zmr)) .lt. abs(Im(zm1r)) ) then zmr = ToComplex(Re(zmr),-Im(zm1r)) else zm1r = ToComplex(Re(zm1r),-Im(zmr)) endif if ( abs(Im(zpr)) .lt. abs(Im(zp1r)) ) then zpr = ToComplex(Re(zpr),-Im(zp1r)) else zp1r = ToComplex(Re(zp1r),-Im(zpr)) endif endif else zmr = zm zm1r = zm1 zpr = zp zp1r = zp1 endif call ffc1lg(cs1,zm,zm1,zmr,zm1r,-1,ier) call ffc1lg(cs2,zp,zp1,zpr,zp1r,+1,ier) cb0p = -clogmm/2 + cs1 + cs2 smax = max(absc(clogmm)/2,absc(cs1),absc(cs2)) if ( absc(cb0p) .lt. xloss*smax ) then call ffwarn(7,ier,absc(cb0p),smax) endif goto 990 * -#] normal case: * #] unequal nonzero masses: * #[ debug: 990 continue * #] debug: *###] ffcb0p: end *###[ ffc1lg: subroutine ffc1lg(cs,z,z1,zr,z1r,is,ier) ***#[*comment:*********************************************************** * * * Calculate the potentially unstable combination -1-z*log(1-1/z) * * = sum_{n=1} 1/(n+1) z^{-n}. * * * * Input z,z1 complex root, z1=1-z * * zr,z1r complex root for Im(p^2)=0, z1r=1-zr * * is integer -1: roots are z-, +1: z+ * * * * Output cs complex see above * * ier integer usual error flag * * * ***#]*comment:*********************************************************** * #[ declarations: implicit none * * arguments * integer is,ier ComplexType cs,z,z1,zr,z1r * * local variables * RealType absc ComplexType c,zfflog * * common blocks * #include "ff.h" * * statement function * absc(c) = abs(Re(c)) + abs(Im(c)) * * #] declarations: * #[ work: if ( 1 .lt. xclogm*absc(z) ) then cs = 0 elseif ( 1 .lt. precc*absc(z) ) then cs = 1/(2*z) elseif ( 1 .gt. 2*xloss*absc(z) ) then * * normal case * cs = -1 - z*zfflog(-z1/z,0,czero,ier) * * check analytical continuation for Im(p^2) -> 0 * if ( z.ne.zr .or. z1.ne.z1r ) then c = -z1r/zr if ( Re(c).lt.0 ) then * check whetehr we chose the correct continuation if ( (Im(c).gt.0 .or. Im(c).eq.0 .and. + is.eq.+1) .and. Im(-z1/z).lt.0 ) then cs = cs - c2ipi*z elseif ( (Im(c).lt.0 .or. Im(c).eq.0 .and. + is.eq.-1) .and. Im(-z1/z).gt.0 ) then cs = cs + c2ipi*z endif endif endif if ( absc(cs) .lt. xloss ) call ffwarn(8,ier,absc(cs),1D0) else * * Taylor expansion * call ffcayl(cs,1/z,xninv(2),29,ier) endif * #] work: *###] ffc1lg: end *###[ ffcot2: subroutine ffcot2(cpiDpj,cp,cma,cmb,cmap,cmbp,cmamb,ier) ***#[*comment:*********************************************************** * * * Store the 3 dotproducts in the common block ffdot. * * * * Input: see ffxc0p * * * ***#]*comment:*********************************************************** * #[ declarations: implicit none * * arguments * integer ier ComplexType cpiDpj(3,3),cp,cma,cmb,cmap,cmbp,cmamb * * local variables * integer ier1 RealType absc,xmax ComplexType c * * common blocks * #include "ff.h" * * statement function * absc(c) = abs(Re(c)) + abs(Im(c)) * #] declarations: * #[ work: ier1 = ier cpiDpj(1,1) = cma cpiDpj(2,2) = cmb cpiDpj(3,3) = cp if ( absc(cmap) .lt. absc(cmbp) ) then cpiDpj(1,2) = (cmap + cmb)/2 else cpiDpj(1,2) = (cmbp + cma)/2 endif cpiDpj(2,1) = cpiDpj(1,2) xmax = min(absc(cma),absc(cmb))/2 if ( absc(cmamb) .lt. absc(cmbp) ) then cpiDpj(1,3) = (-cmamb - cp)/2 else cpiDpj(1,3) = (cmbp - cma)/2 endif cpiDpj(3,1) = cpiDpj(1,3) xmax = min(absc(cma),absc(cp))/2 if ( absc(cmamb) .lt. absc(cmap) ) then cpiDpj(2,3) = (-cmamb + cp)/2 else cpiDpj(2,3) = (-cmap + cmb)/2 endif cpiDpj(3,2) = cpiDpj(2,3) xmax = min(absc(cmb),absc(cp))/2 ier = ier1 * #] work: *###] ffcot2: end diff --git a/Looptools/B/ffxb0.F b/Looptools/B/ffxb0.F --- a/Looptools/B/ffxb0.F +++ b/Looptools/B/ffxb0.F @@ -1,968 +1,969 @@ #include "externals.h" #include "types.h" *###[ ffxb0: subroutine ffxb0(cb0,xp,xma,xmb,ier) ***#[*comment:*********************************************************** * * * Calculates the the two-point function (cf 't Hooft and Veltman) * * we include an overall factor 1/(i*pi^2) relative to FormF * * * * Input: xp (real) k2, in B&D metric * * xma (real) mass2 * * xmb (real) mass2 * * * * Output: cb0 (complex) B0, the two-point function, * * ier (integer) # of digits lost, if >=100: error * * * * Calls: ffxb0p * * * ***#]*comment:*********************************************************** * #[ declarations: implicit none * * arguments * integer ier ComplexType cb0 RealType xp,xma,xmb * * local variables * ComplexType cb0p RealType dmamb,dmap,dmbp,xm * * common blocks * #include "ff.h" * * #] declarations: * #[ get differences: dmamb = xma - xmb dmap = xma - xp dmbp = xmb - xp * #] get differences: * #[ calculations: call ffxb0p(cb0p,xp,xma,xmb,dmap,dmbp,dmamb,ier) if ( xma .eq. 0 ) then if ( xmb .eq. 0 ) then xm = 1D0 else xm = xmb**2 endif elseif ( xmb .eq. 0 ) then xm = xma**2 else xm = xma*xmb endif if ( mudim .ne. 0 ) xm = xm/mudim**2 if ( abs(xm) .gt. xalogm ) then cb0 = Re(delta - log(xm)/2D0) - cb0p else call fferr(4,ier) cb0 = Re(delta) - cb0p endif * #] calculations: *###] ffxb0: end *###[ ffxb0p: subroutine ffxb0p(cb0p,xp,xma,xmb,dmap,dmbp,dmamb,ier) ***#[*comment:*********************************************************** * * * calculates the two-point function (see 't Hooft and * * Veltman) for all possible cases: masses equal, unequal, * * equal to zero. * * * * Input: xp (real) p.p, in B&D metric * * xma (real) mass2, * * xmb (real) mass2, * * dm[ab]p (real) xm[ab] - xp * * dmamb (real) xma - xmb * * * * Output: cb0p (complex) B0, the two-point function, minus * * log(xm1*xm2)/2, delta and ipi^2 * * ier (integer) 0=ok, 1=numerical problems, 2=error * * * * Calls: ffxb0q. * * * ***#]*comment:*********************************************************** * #[ declarations: implicit none * * arguments * integer ier ComplexType cb0p RealType xp,xma,xmb,dmap,dmbp,dmamb * * local variables * integer i,initeq,initn1,jsign RealType ax,ay,ffbnd, + xprceq,bdeq01,bdeq05,bdeq11,bdeq17, + xprcn1,bdn101,bdn105,bdn110,bdn115, + xprnn2,bdn205,bdn210,bdn215,bdn220, + xprcn3,bdn301,bdn305,bdn310,bdn315, + xprcn5,bdn501,bdn505,bdn510,bdn515, + absc RealType xm,dmp,xm1,xm2,dm1m2,dm1p, + dm2p,s,s1,s1a,s1b,s1p,s2,s2a,s2b,s2p,x,y,som, + xlam,slam,xlogmm,alpha,alph1,xnoe,xpneq(30), + xpnn1(30),xx,xtel,dfflo1 ComplexType cs2a,cs2b,cs2p,c,cx external ffbnd,dfflo1 save initeq,initn1,xpneq,xpnn1, + xprceq,bdeq01,bdeq05,bdeq11,bdeq17, + xprcn1,bdn101,bdn105,bdn110,bdn115, + xprnn2,bdn205,bdn210,bdn215,bdn220, + xprcn3,bdn301,bdn305,bdn310,bdn315, + xprcn5,bdn501,bdn505,bdn510,bdn515 * * common blocks * #include "ff.h" * * data * data xprceq /-1D0/ data xprcn1 /-1D0/ data xprnn2 /-1D0/ data xprcn3 /-1D0/ data xprcn5 /-1D0/ data initeq /0/ data initn1 /0/ * * statement function * absc(c) = abs(Re(c)) + abs(Im(c)) * #] declarations: * #[ fill some dotproducts: if ( ldot ) then call ffdot2(fpij2,xp,xma,xmb,dmap,dmbp,dmamb,ier) endif * #] fill some dotproducts: * #[ which case: * * sort according to the type of masscombination encountered: * 100: both masses zero, 200: one equal to zero, 300: both equal * 400: rest. * if ( xma .eq. 0 ) then if ( xmb .eq. 0 ) then goto 100 endif xm = xmb dmp = dmbp goto 200 endif if ( xmb .eq. 0 ) then xm = xma dmp = dmap goto 200 elseif ( dmamb .eq. 0 ) then xm = xma dmp = dmap goto 300 elseif ( xma .gt. xmb ) then xm2 = xma xm1 = xmb dm1m2 = -dmamb dm1p = dmbp dm2p = dmap else xm1 = xma xm2 = xmb dm1m2 = dmamb dm1p = dmap dm2p = dmbp endif goto 400 * #] which case: * #[ both masses equal to zero: 100 continue if ( xp .lt. -xalogm ) then cb0p = log(-xp) - 2 elseif ( xp .gt. xalogm ) then cb0p = ToComplex( Re(log(xp) - 2), Re(-pi) ) else cb0p = 0 call fferr(7,ier) endif return * #] both masses equal to zero: * #[ one mass equal to zero: 200 continue * * special case xp = 0 * if ( xp .eq. 0 ) then cb0p = -1 goto 990 * * special case xp = xm * elseif ( dmp.eq.0 ) then cb0p = -2 goto 990 endif * * Normal case: * s1 = xp/xm if ( abs(s1) .lt. xloss ) then s = dfflo1(s1,ier) else s = log(abs(dmp/xm)) endif s = -s*dmp/xp cb0p = s - 2 if ( xp .gt. xm ) + cb0p = cb0p - ToComplex(0D0,-(dmp/xp)*pi) goto 990 * #] one mass equal to zero: * #[ both masses equal: 300 continue * * Both masses are equal. Not only this speeds up things, some * cancellations have to be avoided as well. * * first a special case * if ( abs(xp) .lt. 8*xloss*xm ) then * -#[ taylor expansion: * * a Taylor expansion seems appropriate as the result will go * as k^2 but seems to go as 1/k !! * *--#[ data and bounds: if ( initeq .eq. 0 ) then initeq = 1 xpneq(1) = 1D0/6D0 do 1 i=2,30 xpneq(i) = - xpneq(i-1)*Re(i-1)/Re(2*(2*i+1)) 1 continue endif if (xprceq .ne. precx ) then * * calculate the boundaries for the number of terms to be * included in the taylorexpansion * xprceq = precx bdeq01 = ffbnd(1,1,xpneq) bdeq05 = ffbnd(1,5,xpneq) bdeq11 = ffbnd(1,11,xpneq) bdeq17 = ffbnd(1,17,xpneq) endif *--#] data and bounds: x = -xp/xm ax = abs(x) if ( ax .gt. bdeq17 ) then som = x*(xpneq(18) + x*(xpneq(19) + x*(xpneq(20) + + x*(xpneq(21) + x*(xpneq(22) + x*(xpneq(23) + + x*(xpneq(24) + x*xpneq(25) ))))))) else som = 0 endif if ( ax .gt. bdeq11 ) then som = x*(xpneq(12) + x*(xpneq(13) + x*(xpneq(14) + + x*(xpneq(15) + x*(xpneq(16) + x*(xpneq(17) + som )))) + )) endif if ( ax .gt. bdeq05 ) then som = x*(xpneq(6) + x*(xpneq(7) + x*(xpneq(8) + x*( + xpneq(9) + x*(xpneq(10) + x*(xpneq(11) + som )))))) endif if ( ax .gt. bdeq01 ) then som = x*(xpneq(2) + x*(xpneq(3) + x*(xpneq(4) + x*( + xpneq(5) + som )))) endif cb0p = x*(xpneq(1)+som) goto 990 * -#] taylor expansion: endif * -#[ normal case: * * normal case * call ffxlmb(xlam,-xp,-xm,-xm,dmp,dmp,0D0) if ( xlam .ge. 0 ) then * cases 1,2 and 4 slam = sqrt(xlam) s2a = dmp + xm s2 = s2a + slam if ( abs(s2) .gt. xloss*slam ) then * looks fine jsign = 1 else s2 = s2a - slam jsign = -1 endif ax = abs(s2/(2*xm)) if ( ax .lt. xalogm ) then s = 0 elseif( ax-1 .lt. .1 .and. s2 .gt. 0 ) then * In this case a quicker and more accurate way is to * calculate log(1-x). s2 = (xp - slam) * the following line is superfluous. s = -slam/xp*dfflo1(s2/(2*xm),ier) else * finally the normal case s = -slam/xp*log(ax) if ( jsign .eq. -1 ) s = -s endif if ( xp .gt. 2*xm ) then * in this case ( xlam>0, so xp>(2*m)^2) ) there also * is an imaginary part y = -pi*slam/xp else y = 0 endif else * the root is complex (k^2 between 0 and (2*m1)^2) slam = sqrt(-xlam) s = 2*slam/xp*atan2(xp,slam) y = 0 endif xx = s - 2 cb0p = ToComplex(Re(xx),Re(y)) goto 990 * -#] normal case: * * #] both masses equal: * #[ unequal nonzero masses: * -#[ get log(xm2/xm1): 400 continue x = xm2/xm1 if ( 1 .lt. xalogm*x ) then call fferr(8,ier) xlogmm = 0 elseif ( abs(x-1) .lt. xloss ) then xlogmm = dfflo1(dm1m2/xm1,ier) else xlogmm = log(x) endif * -#] get log(xm2/xm1): * -#[ xp = 0: * * first a special case * if ( xp .eq. 0 ) then s2 = ((xm2+xm1) / dm1m2)*xlogmm s = - s2 - 2 * save the factor 1/2 for the end if ( abs(s) .lt. xloss*2 ) then * Taylor expansions: choose which one x = dm1m2/xm1 ax = abs(x) if ( ax .lt. .15 .or. precx .gt. 1.E-8 .and. ax + .lt. .3 ) then * * This is the simple Taylor expansion 'n1' * *--#[ data and bounds: * get the coefficients of the taylor expansion if ( initn1 .eq. 0 ) then initn1 = 1 do 410 i = 1,30 - 410 xpnn1(i) = Re(i)/Re((i+1)*(i+2)) + xpnn1(i) = Re(i)/Re((i+1)*(i+2)) + 410 continue endif * determine the boundaries for 1,5,10,15 terms if ( xprcn1 .ne. precx ) then xprcn1 = precx bdn101 = ffbnd(1,1,xpnn1) bdn105 = ffbnd(1,5,xpnn1) bdn110 = ffbnd(1,10,xpnn1) bdn115 = ffbnd(1,15,xpnn1) endif *--#] data and bounds: * calculate: if ( ax .gt. bdn115 ) then s = x*(xpnn1(16) + x*(xpnn1(17) + x*(xpnn1(18) + + x*(xpnn1(19) + x*xpnn1(20) )))) else s = 0 endif if ( ax .gt. bdn110 ) then s = x*(xpnn1(11) + x*(xpnn1(12) + x*(xpnn1(13) + + x*(xpnn1(14) + x*xpnn1(15) + s)))) endif if ( ax .gt. bdn105 ) then s = x*(xpnn1(6) + x*(xpnn1(7) + x*(xpnn1(8) + x* + (xpnn1(9) + x*(xpnn1(10) + s))))) endif if ( ax .gt. bdn101 ) then s = x*(xpnn1(2) + x*(xpnn1(3) + x*(xpnn1(4) + x* + (xpnn1(5) +s)))) endif s = x*x*(xpnn1(1) + s) else * * This is the more complicated Taylor expansion 'fc' * * #[ bounds: * determine the boundaries for 1,5,10,15 terms for * the exponential taylor expansion, assuming it * starts at n=2. * if ( xprnn2 .ne. precx ) then xprnn2 = precx bdn205 = ffbnd(4,5,xinfac) bdn210 = ffbnd(4,10,xinfac) bdn215 = ffbnd(4,15,xinfac) bdn220 = ffbnd(4,20,xinfac) endif * #] bounds: * calculate: y = 2*x/(2-x) ay = abs(y) if ( ay .gt. bdn220 ) then s = y*(xinfac(19) + y*(xinfac(20) + y*(xinfac( + 21) + y*(xinfac(22) + y*xinfac( + 23) )))) else s = 0 endif if ( ay .gt. bdn215 ) then s = y*(xinfac(14) + y*(xinfac(15) + y*(xinfac( + 16) + y*(xinfac(17) + y*(xinfac( + 18) + s))))) endif if ( ay .gt. bdn210 ) then s = y*(xinfac(9) + y*(xinfac(10) + y*(xinfac(11) + + y*(xinfac(12) + y*(xinfac(13) + s))))) endif if ( ay .gt. bdn205 ) then s = y*(xinfac(5) + y*(xinfac(6) + y*(xinfac(7) + + y*(xinfac(8) + s)))) endif s = (1-x)*y**4*(xinfac(4)+s) s = x*y**2*(1+y)/12 - s s = - 2*dfflo1(s,ier)/y endif endif cb0p = s/2 goto 990 endif * -#] xp = 0: * -#[ normal case: * * proceeding with the normal case * call ffxlmb(xlam,-xp,-xm2,-xm1,dm2p,dm1p,dm1m2) if ( xlam .gt. 0 ) then * cases k^2 < -(m2+m1)^2 or k^2 > -(m2-m1)^2: *--#[ first try: * first try the normal way slam = sqrt(xlam) s2a = dm2p + xm1 s2 = s2a + slam if ( abs(s2) .gt. xloss*slam ) then * looks fine jsign = 1 else s2 = s2a - slam jsign = -1 endif s2 = s2**2/(4*xm1*xm2) if ( abs(s2) .lt. xalogm ) then call fferr(9,ier) s2 = 0 elseif ( abs(s2-1) .lt. xloss ) then if ( jsign.eq.1 ) then s2 = -slam*(s2a+slam)/(2*xm1*xm2) s2 = -slam/(2*xp)*dfflo1(s2,ier) else s2 = +slam*(s2a-slam)/(2*xm1*xm2) s2 = +slam/(2*xp)*dfflo1(s2,ier) endif else s2 = -slam/(2*xp)*log(s2) if ( jsign .eq. -1 ) s2 = -s2 endif s1 = -dm1m2*xlogmm/(2*xp) xx = s1+s2-2 *--#] first try: if ( abs(xx) .lt. xloss*max(abs(s1),abs(s2)) ) then *--#[ second try: * this is unacceptable, try a better solution s1a = dm1m2 + slam if ( abs(s1a) .gt. xloss*slam ) then * (strangely) this works s1 = -s1a/(2*xp) else * by division a more accurate form can be found s1 = ( -xp/2 + xm1 + xm2 ) / ( slam - dm1m2 ) endif s1 = s1*xlogmm if ( abs(xp) .lt. xm2 ) then s2a = xp - dm1m2 else s2a = xm2 - dm1p endif s2 = s2a - slam if ( abs(s2) .gt. xloss*slam ) then * at least reasonable s2 = s2 / (2*xm2) else * division again s2 = (2*xp) / (s2a+slam) endif if ( abs(s2) .lt. .1 ) then * choose a quick way to get the logarithm s2 = dfflo1(s2,ier) else s2a = abs(1-s2) s2 = log(s2a) endif s2 = -(slam/xp)*s2 xx = s1 + s2 - 2 *--#] second try: if ( abs(xx) .lt. xloss**2*max(abs(s1),abs(s2)) ) then *--#[ third try: * (we accept two times xloss because that's the same * as in this try) * A Taylor expansion might work. We expand * inside the logs. Only do the necessary work. * alpha = slam/(slam-dm1m2) alph1 = -dm1m2/(slam-dm1m2) * * First s1: * s1p = s1 - 2*alph1 if ( abs(s1p) .lt. xloss*abs(s1) ) then * -#[ bounds: * determine the boundaries for 1,5,10,15 terms if ( xprcn3 .ne. precx ) then xprcn3 = precx bdn301 = ffbnd(3,1,xinfac) bdn305 = ffbnd(3,5,xinfac) bdn310 = ffbnd(3,10,xinfac) bdn315 = ffbnd(3,15,xinfac) endif * -#] bounds: xnoe = -xp + 2*xm1 + 2*xm2 x = 4*dm1m2/xnoe ax = abs(x) if ( ax .gt. bdn310 ) then s1a = x*(xinfac(13) + x*(xinfac(14) + x*( + xinfac(15) + x*(xinfac(16) + x* + xinfac(17) )))) else s1a = 0 endif if ( ax .gt. bdn305 ) then s1a = x*(xinfac(8) + x*(xinfac(9) + x*( + xinfac(10) + x*(xinfac(11) + x*( + xinfac(12) + s1a))))) endif if ( ax .gt. bdn301 ) then s1a = x*(xinfac(4) + x*(xinfac(5) + x*( + xinfac(6) + x*(xinfac(7) + s1a)))) endif s1a = x**3 *(xinfac(3) + s1a) *xm2/xm1 s1b = dm1m2*(4*dm1m2**2 - xp*(4*xm1-xp))/ + (xm1*xnoe**2) s1p = s1b - s1a s1p = xnoe*dfflo1(s1p,ier)/(slam - dm1m2)/2 endif * * next s2: * s2p = s2 - 2*alpha if ( abs(s2p) .lt. xloss*abs(s2) ) then * -#[ bounds: * determine the boundaries for 1,5,10,15 terms if ( xprcn5 .ne. precx ) then xprcn5 = precx bdn501 = ffbnd(4,1,xinfac) bdn505 = ffbnd(4,5,xinfac) bdn510 = ffbnd(4,10,xinfac) bdn515 = ffbnd(4,15,xinfac) endif * -#] bounds: xnoe = slam - dm1m2 x = 2*xp/xnoe ax = abs(x) * do not do the Taylor expansion if ( ax .gt. bdn515 ) goto 495 if ( ax .gt. bdn510 ) then s2a = x*(xinfac(14) + x*(xinfac(15) + x*( + xinfac(16) + x*(xinfac(17) + x* + xinfac(18) )))) else s2a = 0 endif if ( ax .gt. bdn505 ) then s2a = x*(xinfac(9) + x*(xinfac(10) + x*( + xinfac(11) + x*(xinfac(12) + x*( + xinfac(13) + s2a))))) endif if ( ax .gt. bdn501 ) then s2a = x*(xinfac(5) + x*(xinfac(6) + x*( + xinfac(7) + x*(xinfac(8) + s2a)))) endif s2a = x**4*(xinfac(4)+s2a)*(1-2*xp/(xnoe+xp)) s2b = -2*xp**3 *(-2*xp - xnoe)/(3*(xnoe+xp)* + xnoe**3) s2p = s2b - s2a s2p = -slam/xp*dfflo1(s2p,ier) endif * * finally ... * 495 xx = s1p + s2p *--#] third try: endif endif if ( xp .gt. xm1+xm2 ) then *--#[ imaginary part: * in this case ( xlam>0, so xp>(m1+m2)^2) ) there also * is an imaginary part y = -pi*slam/xp else y = 0 *--#] imaginary part: endif else * the root is complex (k^2 between -(m1+m2)^2 and -(m2-m1)^2) *--#[ first try: slam = sqrt(-xlam) xnoe = dm2p + xm1 s1 = -(dm1m2/(2*xp))*xlogmm s2 = (slam/xp)*atan2(slam,xnoe) xx = s1 + s2 - 2 *--#] first try: * 13 Apr 11: added x .ne. 0 check to safeguard against div by zero x = 2*xp*xnoe if ( x .ne. 0 .and. & abs(xx) .lt. xloss**2*max(abs(s1),abs(s2)) ) then *--#[ second try: * Again two times xloss as we'll accept that in the next * step as well. * xtel = dm1m2**2 - xp**2 alpha = -xlam/x alph1 = xtel/x * * try a taylor expansion on the terms. First s1: * s1p = s1 - 2*alph1 if ( abs(s1p) .lt. xloss*abs(s1) ) then * -#[ bounds: * determine the boundaries for 1,5,10,15 terms if ( xprcn3 .ne. precx ) then xprcn3 = precx bdn301 = ffbnd(3,1,xinfac) bdn305 = ffbnd(3,5,xinfac) bdn310 = ffbnd(3,10,xinfac) bdn315 = ffbnd(3,15,xinfac) endif * -#] bounds: x = 2*xtel/(dm1m2*xnoe) ax = abs(x) * do not do the Taylor expansion if ( ax .gt. bdn315 ) goto 590 if ( ax .gt. bdn310 ) then s1a = x*(xinfac(13) + x*(xinfac(14) + x*( + xinfac(15) + x*(xinfac(16) + x* + xinfac(17) )))) else s1a = 0 endif if ( ax .gt. bdn305 ) then s1a = x*(xinfac(8) + x*(xinfac(9) + x*( + xinfac(10) + x*(xinfac(11) + x*( + xinfac(12) + s1a))))) endif if ( ax .gt. bdn301 ) then s1a = x*(xinfac(4) + x*(xinfac(5) + x*( + xinfac(6) + x*(xinfac(7) + s1a)))) endif s1a = x**3 *(xinfac(3) + s1a) *xm2/xm1 s1b = (dm1m2**3*(dm1m2**2-2*xp*xm1) + xp**2*(4* + dm1m2*xm1**2-dm1m2**2*(dm1m2+2*xm1))-2*xm2* + xp**3*(dm1m2+xp))/(xm1*dm1m2**2*xnoe**2) s1p = s1b - s1a s1p = -dm1m2*dfflo1(s1p,ier)/(2*xp) endif * * next s2: * 590 continue s2p = s2 - 2*alpha if ( abs(s2p) .lt. xloss*abs(s2) ) then * -#[ bounds: * determine the boundaries for 1,5,10,15 terms if ( xprcn3 .ne. precx ) then xprcn3 = precx bdn301 = ffbnd(3,1,xinfac) bdn305 = ffbnd(3,5,xinfac) bdn310 = ffbnd(3,10,xinfac) bdn315 = ffbnd(3,15,xinfac) endif * -#] bounds: cx = ToComplex(0D0,-slam/xnoe) ax = absc(cx) if ( ax .gt. bdn315 ) goto 600 if ( ax .gt. bdn310 ) then cs2a = cx*(Re(xinfac(13)) + cx*(Re(xinfac(14 + )) + cx*(Re(xinfac(15)) + cx*(Re(xinfac(16 + )) + cx*(Re(xinfac(17))))))) else cs2a = 0 endif if ( ax .gt. bdn305 ) then cs2a = cx*(Re(xinfac(8)) + cx*(Re(xinfac(9)) + + cx*(Re(xinfac(10)) + cx*(Re(xinfac(11)) + + cx*(Re(xinfac(12)) + cs2a))))) endif if ( ax .gt. bdn301 ) then cs2a = cx*(Re(xinfac(4)) + cx*(Re(xinfac(5)) + + cx*(Re(xinfac(6)) + cx*(Re(xinfac(7)) + + cs2a)))) endif cs2a = cx**3*(Re(xinfac(3))+cs2a)* + ToComplex(Re(xnoe),Re(slam)) cs2b = ToComplex(Re(xnoe-xlam/xnoe/2), + -Re(slam**3/xnoe**2/2)) cs2p = cs2b + cs2a s2p = slam*atan2(Im(cs2p),Re(cs2p))/xp endif 600 continue xx = s1p + s2p *--#] second try: endif y = 0 endif cb0p = ToComplex(Re(xx),Re(y)) goto 990 * -#] normal case: * #] unequal nonzero masses: * #[ debug: 990 continue * #] debug: *###] ffxb0p: end *###[ ffxlmb: subroutine ffxlmb(xlambd,a1,a2,a3,a12,a13,a23) ***#[*comment:*********************************************************** * calculate in a numerically stable way * * lambda(a1,a2,a3) = * * a1**2 + a2**2 + a3**2 - 2*a2*a3 - 2*a3*a1 - 2*a1*a2 * * aij = ai - aj are required for greater accuracy at times * ***#]*comment:*********************************************************** * #[ declarations: implicit none * * arguments * RealType xlambd,a1,a2,a3,a12,a13,a23 * * local variables * RealType aa1,aa2,aa3,a,aff,asq * * common blocks * #include "ff.h" * #] declarations: * #[ calculations: aa1 = abs(a1) aa2 = abs(a2) aa3 = abs(a3) * * first see if there are input parameters with opposite sign: * if ( a1 .lt. 0 .and. a2 .gt. 0 .or. + a1 .gt. 0 .and. a2 .lt. 0 ) then goto 12 elseif ( a1 .lt. 0 .and. a3 .gt. 0 .or. + a1 .gt. 0 .and. a3 .lt. 0 ) then goto 13 * * all have the same sign, choose the smallest 4*ai*aj term * elseif ( aa1 .gt. aa2 .and. aa1 .gt. aa3 ) then goto 23 elseif ( aa2 .gt. aa3 ) then goto 13 else goto 12 endif 12 continue if ( aa1 .gt. aa2 ) then a = a13 + a2 else a = a1 + a23 endif aff = 4*a1*a2 goto 100 13 continue if ( aa1 .gt. aa3 ) then a = a12 + a3 else a = a1 - a23 endif aff = 4*a1*a3 goto 100 23 continue if ( aa2 .gt. aa3 ) then a = a12 - a3 else a = a13 - a2 endif aff = 4*a2*a3 100 continue asq = a**2 xlambd = asq - aff * #] calculations: *###] ffxlmb: end *###[ ffclmb: subroutine ffclmb(clambd,cc1,cc2,cc3,cc12,cc13,cc23) ***#[*comment:*********************************************************** * calculate in cc numerically stable way * * lambda(cc1,cc2,cc3) = * * cc1**2 + cc2**2 + cc3**2 - 2*cc2*cc3 - 2*cc3*cc1 - 2*cc1*cc2 * * cij = ci - cj are required for greater accuracy at times * * ier is the usual error flag. * ***#]*comment:*********************************************************** * #[ declarations: implicit none * * arguments * ComplexType clambd,cc1,cc2,cc3,cc12,cc13,cc23 * * local variables * RealType aa1,aa2,aa3,absc ComplexType cc,cff,csq,c * * common blocks * #include "ff.h" * * statement function * absc(c) = abs(Re(c)) + abs(Im(c)) * #] declarations: * #[ calculations (rather old style ...): aa1 = absc(cc1) aa2 = absc(cc2) aa3 = absc(cc3) * * first see if there are input parameters with opposite sign: * if ( Re(cc1) .lt. 0 .and. Re(cc2) .gt. 0 .or. + Re(cc1) .gt. 0 .and. Re(cc2) .lt. 0 ) then goto 12 elseif ( Re(cc1) .lt. 0 .and. Re(cc3) .gt. 0 .or. + Re(cc1) .gt. 0 .and. Re(cc3) .lt. 0 ) then goto 13 * * all have the same sign, choose the smallest 4*ci*cj term * elseif ( aa1 .gt. aa2 .and. aa1 .gt. aa3 ) then goto 23 elseif ( aa2 .gt. aa3 ) then goto 13 else goto 12 endif 12 continue if ( aa1 .gt. aa2 ) then cc = cc13 + cc2 else cc = cc1 + cc23 endif cff = 4*cc1*cc2 goto 100 13 continue if ( aa1 .gt. aa3 ) then cc = cc12 + cc3 else cc = cc1 - cc23 endif cff = 4*cc1*cc3 goto 100 23 continue if ( aa2 .gt. aa3 ) then cc = cc12 - cc3 else cc = cc13 - cc2 endif cff = 4*cc2*cc3 100 continue csq = cc**2 clambd = csq - cff * #] calculations (rather old style ...): *###] ffclmb: end *###[ ffdot2: subroutine ffdot2(piDpj,xp,xma,xmb,dmap,dmbp,dmamb,ier) ***#[*comment:*********************************************************** * * * Store the 3 dotproducts in the common block ffdot. * * * * Input: see ffxb0p * * * ***#]*comment:*********************************************************** * #[ declarations: implicit none * * arguments * integer ier RealType piDpj(3,3),xp,xma,xmb,dmap,dmbp,dmamb * * local variables * integer ier1 * * common blocks * #include "ff.h" * * statement function * * absc(c) = abs(Re(c)) + abs(Im(c)) * #] declarations: * #[ work: ier1 = ier piDpj(1,1) = xma piDpj(2,2) = xmb piDpj(3,3) = xp if ( abs(dmap) .lt. abs(dmbp) ) then piDpj(1,2) = (dmap + xmb)/2 else piDpj(1,2) = (dmbp + xma)/2 endif piDpj(2,1) = piDpj(1,2) if ( abs(dmamb) .lt. abs(dmbp) ) then piDpj(1,3) = (-dmamb - xp)/2 else piDpj(1,3) = (dmbp - xma)/2 endif piDpj(3,1) = piDpj(1,3) if ( abs(dmamb) .lt. abs(dmap) ) then piDpj(2,3) = (-dmamb + xp)/2 else piDpj(2,3) = (-dmap + xmb)/2 endif piDpj(3,2) = piDpj(2,3) ier = ier1 * #] work: *###] ffdot2: end diff --git a/Shower/Dipole/DipoleShowerHandler.cc b/Shower/Dipole/DipoleShowerHandler.cc --- a/Shower/Dipole/DipoleShowerHandler.cc +++ b/Shower/Dipole/DipoleShowerHandler.cc @@ -1,1384 +1,1384 @@ // -*- C++ -*- // // DipoleShowerHandler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the DipoleShowerHandler class. // #include #include "DipoleShowerHandler.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/RefVector.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" // include theses to have complete types #include "Herwig/PDF/MPIPDF.h" #include "Herwig/PDF/MinBiasPDF.h" #include "Herwig/PDF/HwRemDecayer.h" #include "Herwig/Shower/Dipole/Utility/DipolePartonSplitter.h" #include "Herwig/MatrixElement/Matchbox/Base/MergerBase.h" #include "Herwig/MatrixElement/Matchbox/Base/SubtractedME.h" #include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h" #include using namespace Herwig; bool DipoleShowerHandler::firstWarn = true; DipoleShowerHandler::DipoleShowerHandler() : ShowerHandler(), chainOrderVetoScales(true), nEmissions(0), discardNoEmissions(false), firstMCatNLOEmission(false), thePowhegDecayEmission(true), realignmentScheme(0), verbosity(0), printEvent(0), nTries(0), didRadiate(false), didRealign(false), theRenormalizationScaleFreeze(1.*GeV), theFactorizationScaleFreeze(2.*GeV), theDoCompensate(false), theFreezeGrid(500000), theDetuning(1.0), maxPt(ZERO), muPt(ZERO), theInputColouredOffShellInShower(), theZBoundaries(1) {} DipoleShowerHandler::~DipoleShowerHandler() {} IBPtr DipoleShowerHandler::clone() const { return new_ptr(*this); } IBPtr DipoleShowerHandler::fullclone() const { return new_ptr(*this); } void DipoleShowerHandler::cascade(tPVector ) { throw Exception() << "DipoleShowerHandler: Dipoleshower not implemented as second shower." << "Check your setup or contact Herwig authors." << Exception::runerror; } tPPair DipoleShowerHandler::cascade(tSubProPtr sub, XCombPtr, Energy optHardPt, Energy optCutoff) { useMe(); prepareCascade(sub); resetWeights(); if ( !doFSR() && ! doISR() ) return sub->incoming(); eventRecord().clear(); eventRecord().prepare(sub, dynamic_ptr_cast(lastXCombPtr()), newStep(), pdfs(), ShowerHandler::currentHandler()->generator()->currentEvent()->incoming(), firstInteraction(), offShellPartons()); if ( eventRecord().outgoing().empty() && !doISR() ) return sub->incoming(); if ( !eventRecord().incoming().first->coloured() && !eventRecord().incoming().second->coloured() && !doFSR() ) return sub->incoming(); nTries = 0; while ( true ) { try { didRadiate = false; didRealign = false; if ( eventRecord().truncatedShower() ) { throw Exception() << "Inconsistent hard emission set-up in DipoleShowerHandler::cascade. " << "No truncated shower needed with DipoleShowerHandler. Add " << "'set MEMatching:TruncatedShower No' to input file." << Exception::runerror; } hardScales(lastXCombPtr()->lastShowerScale()); if ( verbosity > 1 ) { generator()->log() << "DipoleShowerHandler starting off:\n"; eventRecord().debugLastEvent(generator()->log()); generator()->log() << flush; } unsigned int nEmitted = 0; if ( firstMCatNLOEmission ) { if ( !eventRecord().isMCatNLOHEvent() ) nEmissions = 1; else nEmissions = 0; } if ( !firstMCatNLOEmission ) { doCascade(nEmitted,optHardPt,optCutoff); if ( discardNoEmissions ) { if ( !didRadiate ) throw Veto(); if ( nEmissions ) if ( nEmissions < nEmitted ) throw Veto(); } } else { if ( nEmissions == 1 ) doCascade(nEmitted,optHardPt,optCutoff); } if ( intrinsicPtGenerator ) { if ( eventRecord().incoming().first->coloured() && eventRecord().incoming().second->coloured() ) { SpinOneLorentzRotation rot = intrinsicPtGenerator->kick(eventRecord().incoming(), eventRecord().intermediates()); eventRecord().transform(rot); } } didRealign = realign(); constituentReshuffle(); // Decay and shower any particles that require decaying while ( !eventRecord().decays().empty() ) { map::const_iterator decayIt = eventRecord().decays().begin(); // find the decay to do, one with greatest width and parent showered while(find(eventRecord().outgoing().begin(),eventRecord().outgoing().end(),decayIt->first)== eventRecord().outgoing().end() && find(eventRecord().hard().begin(),eventRecord().hard().end(),decayIt->first)== eventRecord().hard().end()) ++decayIt; assert(decayIt!=eventRecord().decays().end()); PPtr incoming = decayIt->first; eventRecord().currentDecay(decayIt->second); // Use this to record if an emission actually happens bool powhegEmission = !( nEmissions && nEmitted==nEmissions) ? thePowhegDecayEmission : false; // Decay the particle / sort out its pert proc Energy showerScale = eventRecord().decay(incoming, powhegEmission); // Following the decay, the bool powheg emission is updated // to indicate whether or not an emission occurred if ( powhegEmission ) nEmitted += 1; // Check that there is only one particle incoming to the decay assert(eventRecord().currentDecay()->incoming().size()==1); // Prepare the event record for the showering of the decay bool needToShower = eventRecord().prepareDecay(eventRecord().currentDecay(), offShellPartons()); // Only need to shower if we have coloured outgoing particles if ( needToShower ) { // The decays currently considered produce a maximum of 2 chains (with powheg emission) // so all dipole should have the same scale as returned by the decay function. assert( eventRecord().chains().size() <= 2 ); for ( auto & ch : eventRecord().chains()) { for ( auto & dip : ch.dipoles()) { assert ( showerScale > ZERO ); dip.leftScale( showerScale ); dip.rightScale( showerScale ); } } // Perform the cascade doCascade(nEmitted,optHardPt,optCutoff,true); // Do the constituent mass shell reshuffling decayConstituentReshuffle(eventRecord().currentDecay()); } // Update the decays, adding any decays and updating momenta eventRecord().updateDecays(eventRecord().currentDecay()); eventRecord().decays().erase(decayIt); } break; } catch (RedoShower&) { resetWeights(); if ( ++nTries > maxtry() ) throw ShowerTriesVeto(maxtry()); eventRecord().clear(); eventRecord().prepare(sub, dynamic_ptr_cast(lastXCombPtr()), newStep(), pdfs(), ShowerHandler::currentHandler()->generator()->currentEvent()->incoming(), firstInteraction(), offShellPartons()); continue; } catch (...) { throw; } } tPPair incoming=eventRecord().fillEventRecord(newStep(),firstInteraction(),didRealign); setDidRunCascade(true); return incoming; } // Reshuffle the outgoing partons from the hard process onto their constituent mass shells void DipoleShowerHandler::constituentReshuffle() { if ( constituentReshuffler && ShowerHandler::currentHandler()->retConstituentMasses() ) { if ( eventRecord().decays().empty() ) { constituentReshuffler->reshuffle(eventRecord().outgoing(), eventRecord().incoming(), eventRecord().intermediates()); return; } else { PList decaying; for(auto const & dec : eventRecord().decays()) decaying.push_back(dec.first); constituentReshuffler->hardProcDecayReshuffle( decaying, eventRecord().outgoing(), eventRecord().hard(), eventRecord().incoming(), eventRecord().intermediates()); } } // After reshuffling the hard process, the decays need to be updated // as this is not done in reshuffle vector > decays; for(auto const & dec : eventRecord().decays() ) decays.push_back({dec.first,dec.second}); for(auto const & dec : decays) { PPtr unstable = dec.first; PList::iterator pos = find(eventRecord().intermediates().begin(), eventRecord().intermediates().end(), dec.first); // Update the PPtr in theDecays if(pos!=eventRecord().intermediates().end()) { unstable = *pos; while(!unstable->children().empty()) { unstable = unstable->children()[0]; } eventRecord().decays().erase(dec.first); eventRecord().decays()[unstable] = dec.second; // Update the momenta of any other particles in the decay chain // (for externally provided events) if ( !(eventRecord().decays()[unstable]->outgoing().empty()) ) eventRecord().updateDecayChainMom( unstable , eventRecord().decays()[unstable]); } else { if ( !(eventRecord().decays()[unstable]->outgoing().empty()) ) { // Update the momenta of any other particles in the decay chain // (for externally provided events) // Note this needs to be done for all decaying particles in the // outgoing/hard regardless of whether that particle radiated // or was involved in the reshuffling, this is due to the // transformation performed for IILightKinematics. if ( (find(eventRecord().outgoing().begin(), eventRecord().outgoing().end(), unstable) != eventRecord().outgoing().end()) || (find(eventRecord().hard().begin(), eventRecord().hard().end(), unstable) != eventRecord().hard().end()) ) eventRecord().updateDecayChainMom( unstable , eventRecord().decays()[unstable]); } } } eventRecord().currentDecay(PerturbativeProcessPtr()); } // Reshuffle outgoing partons from a decay process onto their constituent mass shells void DipoleShowerHandler::decayConstituentReshuffle(PerturbativeProcessPtr decayProc) { if ( Debug::level > 2 ){ // Test this function by comparing the // invariant mass of the outgoing decay // systems before and after reshuffling Lorentz5Momentum testOutMomBefore (ZERO,ZERO,ZERO,ZERO); Energy testInvMassBefore = ZERO; for ( auto const & testDecayOutItBefore : decayProc->outgoing() ) { testOutMomBefore += testDecayOutItBefore.first->momentum(); } testInvMassBefore = testOutMomBefore.m(); // decayReshuffle updates both the event record and the decay perturbative process if ( constituentReshuffler && ShowerHandler::currentHandler()->retConstituentMasses()) { constituentReshuffler->decayReshuffle(decayProc, eventRecord().outgoing(), eventRecord().hard(), eventRecord().intermediates()); } Lorentz5Momentum testOutMomAfter (ZERO,ZERO,ZERO,ZERO); Energy testInvMassAfter = ZERO; for ( auto const & testDecayOutItAfter : decayProc->outgoing() ) { testOutMomAfter += testDecayOutItAfter.first->momentum(); } testInvMassAfter = testOutMomAfter.m(); Energy incomingMass = decayProc->incoming()[0].first->momentum().m(); assert( abs(testInvMassBefore-incomingMass)/GeV < 1e-5 ); assert( abs(testInvMassBefore-testInvMassAfter)/GeV < 1e-5); }else{ // decayReshuffle updates both the event record and the decay perturbative process if ( constituentReshuffler && ShowerHandler::currentHandler()->retConstituentMasses() ) { constituentReshuffler->decayReshuffle(decayProc, eventRecord().outgoing(), eventRecord().hard(), eventRecord().intermediates()); } return; } } // Sets the scale of each particle in the dipole chains by finding the smallest //of several upper bound energy scales: the CMEnergy of the event, //the transverse mass of outgoing particles, the hardScale (maxPT or maxQ) //calculated for each dipole (in both configurations) and the veto scale for each particle void DipoleShowerHandler::hardScales(Energy2 muf) { // Initalise maximum pt as max CMEnergy of the event maxPt = generator()->maximumCMEnergy(); if ( restrictPhasespace() ) { // First interaction == hard collision (i.e. not a MPI collision) if ( !hardScaleIsMuF() || !firstInteraction() ) { if ( !eventRecord().outgoing().empty() ) { for ( auto const & p : eventRecord().outgoing() ) maxPt = min(maxPt,p->momentum().mt()); } //Look at any non-coloured outgoing particles in the current subprocess else { assert(!eventRecord().hard().empty()); Lorentz5Momentum phard(ZERO,ZERO,ZERO,ZERO); for ( auto const & p : eventRecord().hard()) phard += p->momentum(); Energy mhard = phard.m(); maxPt = mhard; } maxPt *= hardScaleFactor(); } else { maxPt = hardScaleFactor()*sqrt(muf); } muPt = maxPt; } else { muPt = hardScaleFactor()*sqrt(muf); } for ( auto & ch : eventRecord().chains()) { // Note that minVetoScale is a value for each DipoleChain, not each dipole // It will contain the minimum veto scale from all of the dipoles in the chain Energy minVetoScale = -1.*GeV; for ( auto & dip : ch.dipoles()) { // max scale per config Energy maxFirst = ZERO; Energy maxSecond = ZERO; // Loop over the kernels for the given dipole. // For each dipole configuration, calculate ptMax (or QMax if virtuality ordering) // for each kernel and find the maximum for ( auto const & k : kernels) { pair conf = {true,false}; if ( k->canHandle(dip.index(conf)) ) { // Look in DipoleChainOrdering for this Energy scale = evolutionOrdering()->hardScale(dip.emitter(conf),dip.spectator(conf), dip.emitterX(conf),dip.spectatorX(conf), *k,dip.index(conf)); maxFirst = max(maxFirst,scale); } conf = {false,true}; if ( k->canHandle(dip.index(conf)) ) { Energy scale = evolutionOrdering()->hardScale(dip.emitter(conf),dip.spectator(conf), dip.emitterX(conf),dip.spectatorX(conf), *k,dip.index(conf)); maxSecond = max(maxSecond,scale); } } // Find the maximum value from comparing the maxScale found from maxPt and the vetoScale of the particle if ( dip.leftParticle()->vetoScale() >= ZERO ) { maxFirst = min(maxFirst,sqrt(dip.leftParticle()->vetoScale())); // minVetoScale is a value for each DipoleChain, not each dipole // It contains the minimum veto scale for all the dipoles in the entire DipoleChain if ( minVetoScale >= ZERO ) minVetoScale = min(minVetoScale,sqrt(dip.leftParticle()->vetoScale())); else minVetoScale = sqrt(dip.leftParticle()->vetoScale()); } if ( dip.rightParticle()->vetoScale() >= ZERO ) { maxSecond = min(maxSecond,sqrt(dip.rightParticle()->vetoScale())); if ( minVetoScale >= ZERO ) minVetoScale = min(minVetoScale,sqrt(dip.rightParticle()->vetoScale())); else minVetoScale = sqrt(dip.rightParticle()->vetoScale()); } // Set the emitterScale for both members of each dipole maxFirst = min(maxPt,maxFirst); dip.emitterScale({true,false},maxFirst); maxSecond = min(maxPt,maxSecond); dip.emitterScale({false,true},maxSecond); } // if the smallest veto scale (i.e. from all of the dipoles) // is smaller than the scale calculated for a particular // particle in a particular dipole, // replace the scale with the veto scale if ( !evolutionOrdering()->independentDipoles() && chainOrderVetoScales && minVetoScale >= ZERO ) { for ( auto & dip : ch.dipoles() ) { dip.leftScale(min(dip.leftScale(),minVetoScale)); dip.rightScale(min(dip.rightScale(),minVetoScale)); } } } } Energy DipoleShowerHandler::getWinner(DipoleSplittingInfo& winner, const Dipole& dip, pair conf, Energy optHardPt, Energy optCutoff) { return getWinner(winner,dip.index(conf), dip.emitterX(conf),dip.spectatorX(conf), conf,dip.emitter(conf),dip.spectator(conf), dip.emitterScale(conf),optHardPt,optCutoff); } Energy DipoleShowerHandler::getWinner(SubleadingSplittingInfo& winner, Energy optHardPt, Energy optCutoff) { return getWinner(winner,winner.index(), winner.emitterX(),winner.spectatorX(), winner.configuration(), winner.emitter(),winner.spectator(), winner.startScale(),optHardPt,optCutoff); } Energy DipoleShowerHandler::getWinner(DipoleSplittingInfo& winner, const DipoleIndex& index, double emitterX, double spectatorX, pair conf, tPPtr emitter, tPPtr spectator, Energy startScale, Energy optHardPt, Energy optCutoff) { if ( !index.initialStateEmitter() && !doFSR() ) { winner.didStopEvolving(); return 0.0*GeV; } if ( index.initialStateEmitter() && !doISR() ) { winner.didStopEvolving(); return 0.0*GeV; } // Currently do not split IF dipoles so // don't evaluate them in order to avoid // exceptions in the log if ( index.incomingDecayEmitter() ) { winner.didStopEvolving(); return 0.0*GeV; } DipoleSplittingInfo candidate; candidate.index(index); candidate.configuration(conf); candidate.emitterX(emitterX); candidate.spectatorX(spectatorX); if ( generators().find(candidate.index()) == generators().end() ) getGenerators(candidate.index(),theSplittingReweight); // // NOTE -- needs proper fixing at some point // // For some very strange reason, equal_range gives back // key ranges it hasn't been asked for. This particularly // happens e.g. for FI dipoles of the same kind, but different // PDF (hard vs MPI PDF). I can't see a reason for this, // as DipoleIndex properly implements comparison for equality // and (lexicographic) ordering; for the time being, we // use equal_range, extented by an explicit check for wether // the key is indeed what we wanted. See line after (*) comment // below. // // SW - Update 04/01/2016: Note - This caused a bug for me as I did not // include equality checks on the decay booleans in the == definition pair gens = generators().equal_range(candidate.index()); Energy winnerScale = 0.0*GeV; GeneratorMap::iterator winnerGen = generators().end(); for ( GeneratorMap::iterator gen = gens.first; gen != gens.second; ++gen ) { // (*) see NOTE above if ( !(gen->first == candidate.index()) ) continue; if ( startScale <= gen->second->splittingKinematics()->IRCutoff() ) continue; Energy dScale = gen->second->splittingKinematics()->dipoleScale(emitter->momentum(), spectator->momentum()); // in very exceptional cases happening in DIS if ( std::isnan( double(dScale/MeV) ) ) throw RedoShower(); candidate.scale(dScale); // Calculate the mass of the recoil system // for decay dipoles if ( candidate.index().incomingDecaySpectator() || candidate.index().incomingDecayEmitter() ) { Energy recoilMass = gen->second->splittingKinematics()->recoilMassKin(emitter->momentum(), spectator->momentum()); candidate.recoilMass(recoilMass); } // Store emitter and spectator masses, needed in kinematics if ( candidate.index().emitterData()->mass() != ZERO ) { if ( !candidate.index().offShellEmitter() ) candidate.emitterMass( emitter->nominalMass() ); else candidate.emitterMass( emitter->mass() ); } if ( candidate.index().spectatorData()->mass() != ZERO ) { if ( !candidate.index().offShellSpectator() ) candidate.spectatorMass( spectator->nominalMass() ); else candidate.spectatorMass( spectator->mass() ); } candidate.continuesEvolving(); Energy hardScale = evolutionOrdering()->maxPt(startScale,candidate,*(gen->second->splittingKernel())); Energy maxPossible = gen->second->splittingKinematics()->ptMax(candidate.scale(), candidate.emitterX(), candidate.spectatorX(), candidate, *gen->second->splittingKernel()); Energy ircutoff = optCutoff < gen->second->splittingKinematics()->IRCutoff() ? gen->second->splittingKinematics()->IRCutoff() : optCutoff; if ( maxPossible <= ircutoff ) { continue; } if ( maxPossible >= hardScale ){ candidate.hardPt(hardScale); } else { hardScale = maxPossible; candidate.hardPt(maxPossible); } gen->second->generate(candidate,currentWeights(),optHardPt,optCutoff); Energy nextScale = evolutionOrdering()->evolutionScale( gen->second->lastSplitting(),*(gen->second->splittingKernel())); if ( nextScale > winnerScale ) { winner.fill(candidate); gen->second->completeSplitting(winner); winnerGen = gen; winnerScale = nextScale; } reweight(reweight() * gen->second->splittingWeight()); } if ( winnerGen == generators().end() ) { winner.didStopEvolving(); return 0.0*GeV; } if ( winner.stoppedEvolving() ) return 0.0*GeV; return winnerScale; } void DipoleShowerHandler::doCascade(unsigned int& emDone, Energy optHardPt, Energy optCutoff, const bool decay) { if ( nEmissions ) if ( emDone == nEmissions ) return; DipoleSplittingInfo winner; DipoleSplittingInfo dipoleWinner; while ( eventRecord().haveChain() ) { // allow the dipole chain to be rearranged according to arXiv:1801.06113 - if( _rearrange && ( _rearrangeNEmissions < 0 || _rearrangeNEmissions >= emDone ) ){ + if( _rearrange && ( _rearrangeNEmissions < 0 || _rearrangeNEmissions >= int(emDone) ) ){ eventRecord().currentChain().rearrange(_dipmax,_diplong); } if ( verbosity > 2 ) { generator()->log() << "DipoleShowerHandler selecting splittings for the chain:\n" << eventRecord().currentChain() << flush; } list::iterator winnerDip = eventRecord().currentChain().dipoles().end(); Energy winnerScale = 0.0*GeV; Energy nextLeftScale = 0.0*GeV; Energy nextRightScale = 0.0*GeV; for ( list::iterator dip = eventRecord().currentChain().dipoles().begin(); dip != eventRecord().currentChain().dipoles().end(); ++dip ) { nextLeftScale = getWinner(dipoleWinner,*dip,{true,false},optHardPt,optCutoff); if ( nextLeftScale > winnerScale ) { winnerScale = nextLeftScale; winner = dipoleWinner; winnerDip = dip; } nextRightScale = getWinner(dipoleWinner,*dip,{false,true},optHardPt,optCutoff); if ( nextRightScale > winnerScale ) { winnerScale = nextRightScale; winner = dipoleWinner; winnerDip = dip; } if ( evolutionOrdering()->independentDipoles() ) { Energy dipScale = max(nextLeftScale,nextRightScale); if ( dip->leftScale() > dipScale ) dip->leftScale(dipScale); if ( dip->rightScale() > dipScale ) dip->rightScale(dipScale); } } if ( verbosity > 1 ) { if ( winnerDip != eventRecord().currentChain().dipoles().end() ) generator()->log() << "DipoleShowerHandler selected the splitting:\n" << winner << " for the dipole\n" << (*winnerDip) << flush; else generator()->log() << "DipoleShowerHandler could not select a splitting above the IR cutoff\n" << flush; } // pop the chain if no dipole did radiate if ( winnerDip == eventRecord().currentChain().dipoles().end() ) { eventRecord().popChain(); if ( theEventReweight && eventRecord().chains().empty() ) if ( (theEventReweight->firstInteraction() && firstInteraction()) || (theEventReweight->secondaryInteractions() && !firstInteraction()) ) { double w = theEventReweight->weightCascade(eventRecord().incoming(), eventRecord().outgoing(), eventRecord().hard(),theGlobalAlphaS); reweight(reweight()*w); } continue; } // otherwise perform the splitting // but first see if the emission would produce a configuration in the ME region. if ( theMergingHelper && eventHandler()->currentCollision() && !decay && firstInteraction() ) { if (theMergingHelper->maxLegs()>eventRecord().outgoing().size()+ eventRecord().hard().size() +2){//incoming if (theMergingHelper->mergingScale()emissionProbability() < UseRandom::rnd()) { theMergingHelper->setEmissionProbability(0.); const bool transparent=true; if (transparent) { pair::iterator,list::iterator> tmpchildren; DipoleSplittingInfo tmpwinner=winner; DipoleChain* tmpfirstChain = nullptr; DipoleChain* tmpsecondChain = nullptr; auto New=eventRecord().tmpsplit(winnerDip,tmpwinner, tmpchildren,tmpfirstChain, tmpsecondChain); if (theMergingHelper->matrixElementRegion(New.first, New.second, winnerScale, theMergingHelper->mergingScale())) { optHardPt=winnerScale; continue; } }else{ optHardPt=winnerScale; continue; } } } } if(theMergingHelper&&firstInteraction()) optHardPt=ZERO; didRadiate = true; eventRecord().isMCatNLOSEvent(false); eventRecord().isMCatNLOHEvent(false); pair::iterator,list::iterator> children; DipoleChain* firstChain = nullptr; DipoleChain* secondChain = nullptr; // Note: the dipoles are updated in eventRecord().split(....) after the splitting, // hence the entire cascade is handled in doCascade // The dipole scales are updated in dip->split(....) if ( decay ) winner.isDecayProc( true ); eventRecord().split(winnerDip,winner,children,firstChain,secondChain); assert(firstChain && secondChain); evolutionOrdering()->setEvolutionScale(winnerScale,winner,*firstChain,children); if ( !secondChain->dipoles().empty() ) evolutionOrdering()->setEvolutionScale(winnerScale,winner,*secondChain,children); if ( verbosity > 1 ) { generator()->log() << "DipoleShowerHandler did split the last selected dipole into:\n" << (*children.first) << (*children.second) << flush; } if ( verbosity > 2 ) { generator()->log() << "After splitting the last selected dipole, " << "DipoleShowerHandler encountered the following chains:\n" << (*firstChain) << (*secondChain) << flush; } if ( theEventReweight ) if ( (theEventReweight->firstInteraction() && firstInteraction()) || (theEventReweight->secondaryInteractions() && !firstInteraction()) ) { double w = theEventReweight->weight(eventRecord().incoming(), eventRecord().outgoing(), eventRecord().hard(),theGlobalAlphaS); reweight(reweight()*w); } if ( nEmissions ) if ( ++emDone == nEmissions ) return; } } bool DipoleShowerHandler::realign() { if ( !didRadiate && !intrinsicPtGenerator ) return false; if ( eventRecord().incoming().first->coloured() || eventRecord().incoming().second->coloured() ) { if ( eventRecord().incoming().first->momentum().perp2()/GeV2 < 1e-10 && eventRecord().incoming().second->momentum().perp2()/GeV2 < 1e-10 ) return false; pair inMomenta (eventRecord().incoming().first->momentum(), eventRecord().incoming().second->momentum()); SpinOneLorentzRotation transform((inMomenta.first+inMomenta.second).findBoostToCM()); Axis dir = (transform * inMomenta.first).vect().unit(); Axis rot (-dir.y(),dir.x(),0); double theta = dir.theta(); if ( lastParticles().first->momentum().z() < ZERO ) theta = -theta; transform.rotate(-theta,rot); inMomenta.first = transform*inMomenta.first; inMomenta.second = transform*inMomenta.second; assert(inMomenta.first.z() > ZERO && inMomenta.second.z() < ZERO); Energy2 sHat = (eventRecord().incoming().first->momentum() + eventRecord().incoming().second->momentum()).m2(); pair masses(eventRecord().incoming().first->mass(), eventRecord().incoming().second->mass()); pair qs; if ( !eventRecord().incoming().first->coloured() ) { assert(masses.second == ZERO); qs.first = eventRecord().incoming().first->momentum().z(); qs.second = (sHat-sqr(masses.first))/(2.*(qs.first+sqrt(sqr(masses.first)+sqr(qs.first)))); } else if ( !eventRecord().incoming().second->coloured() ) { assert(masses.first == ZERO); qs.second = eventRecord().incoming().second->momentum().z(); qs.first = (sHat-sqr(masses.second))/(2.*(qs.second+sqrt(sqr(masses.second)+sqr(qs.second)))); } else { assert(masses.first == ZERO && masses.second == ZERO); if ( realignmentScheme == 0 ) { double yX = eventRecord().pX().rapidity(); double yInt = (transform*eventRecord().pX()).rapidity(); double dy = yX-yInt; qs.first = (sqrt(sHat)/2.)*exp(dy); qs.second = (sqrt(sHat)/2.)*exp(-dy); } else if ( realignmentScheme == 1 ) { Energy sS = sqrt((lastParticles().first->momentum() + lastParticles().second->momentum()).m2()); qs.first = eventRecord().fractions().first * sS / 2.; qs.second = eventRecord().fractions().second * sS / 2.; } } double beta = (qs.first-qs.second) / ( sqrt(sqr(masses.first)+sqr(qs.first)) + sqrt(sqr(masses.second)+sqr(qs.second)) ); transform.boostZ(beta); Lorentz5Momentum tmp; if ( eventRecord().incoming().first->coloured() ) { tmp = eventRecord().incoming().first->momentum(); tmp = transform * tmp; eventRecord().incoming().first->set5Momentum(tmp); } if ( eventRecord().incoming().second->coloured() ) { tmp = eventRecord().incoming().second->momentum(); tmp = transform * tmp; eventRecord().incoming().second->set5Momentum(tmp); } eventRecord().transform(transform); return true; } return false; } void DipoleShowerHandler::resetAlphaS(Ptr::tptr as) { for ( auto & k : kernels) { if ( !k->alphaS() ) k->alphaS(as); k->renormalizationScaleFreeze(theRenormalizationScaleFreeze); k->factorizationScaleFreeze(theFactorizationScaleFreeze); } // clear the generators to be rebuild // actually, there shouldn't be any generators // when this happens. generators().clear(); } void DipoleShowerHandler::resetReweight(Ptr::tptr rw) { for ( auto & g : generators() ) g.second->splittingReweight(rw); } void DipoleShowerHandler::getGenerators(const DipoleIndex& ind, Ptr::tptr rw) { bool gotone = false; for ( auto & k : kernels ) { if ( k->canHandle(ind) ) { if ( verbosity > 0 ) { generator()->log() << "DipoleShowerHandler encountered the dipole configuration\n" << ind << " in event number " << eventHandler()->currentEvent()->number() << "\nwhich can be handled by the splitting kernel '" << k->name() << "'.\n" << flush; } gotone = true; Ptr::ptr nGenerator = new_ptr(DipoleSplittingGenerator()); nGenerator->doCompensate(theDoCompensate); nGenerator->splittingKernel(k); if ( renormalizationScaleFactor() != 1. ) nGenerator->splittingKernel()->renormalizationScaleFactor(renormalizationScaleFactor()); if ( factorizationScaleFactor() != 1. ) nGenerator->splittingKernel()->factorizationScaleFactor(factorizationScaleFactor()); if ( !nGenerator->splittingReweight() ) nGenerator->splittingReweight(rw); nGenerator->splittingKernel()->freezeGrid(theFreezeGrid); nGenerator->splittingKernel()->detuning(theDetuning); GeneratorMap::const_iterator equivalent = generators().end(); for ( GeneratorMap::const_iterator eq = generators().begin(); eq != generators().end(); ++eq ) { if ( !eq->second->wrapping() ) if ( k->canHandleEquivalent(ind,*(eq->second->splittingKernel()),eq->first) ) { equivalent = eq; if ( verbosity > 0 ) { generator()->log() << "The dipole configuration " << ind << " can equivalently be handled by the existing\n" << "generator for configuration " << eq->first << " using the kernel '" << eq->second->splittingKernel()->name() << "'\n" << flush; } break; } } if ( equivalent != generators().end() ) { nGenerator->wrap(equivalent->second); } DipoleSplittingInfo dummy; dummy.index(ind); nGenerator->prepare(dummy); generators().insert({ind,nGenerator}); } } if ( !gotone ) { throw Exception() << "DipoleShowerHandler could not " << "find a splitting kernel which is able " << "to handle splittings off the dipole " << ind << ".\n" << "Please check the input files." << Exception::runerror; } } // If needed, insert default implementations of virtual function defined // in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs). void DipoleShowerHandler::doinit() { ShowerHandler::doinit(); if ( theGlobalAlphaS ) resetAlphaS(theGlobalAlphaS); // copy off-shell particle ids before showering from input vector to the // set used in the simulation if ( theColouredOffShellInShower.empty() ) { for(unsigned int ix=0;ixsplittingKinematics()->openZBoundaries() == 0 ) zChoice0 = true; else if ( k->splittingKinematics()->openZBoundaries() == 1 ) zChoice1 = true; else zChoiceOther = true; // either inconsistent or other option which cannot be handled by the matching if ( zChoice0 && zChoice1 ) { zChoiceOther = true; break; } } if ( zChoiceOther ) theZBoundaries = 2; else if ( zChoice1 ) theZBoundaries = 1; else if ( zChoice0 ) theZBoundaries = 0; } void DipoleShowerHandler::dofinish() { ShowerHandler::dofinish(); } void DipoleShowerHandler::doinitrun() { ShowerHandler::doinitrun(); } void DipoleShowerHandler::persistentOutput(PersistentOStream & os) const { os << kernels << theEvolutionOrdering << constituentReshuffler << intrinsicPtGenerator << theGlobalAlphaS << chainOrderVetoScales << nEmissions << discardNoEmissions << firstMCatNLOEmission << thePowhegDecayEmission << realignmentScheme << verbosity << printEvent << ounit(theRenormalizationScaleFreeze,GeV) << ounit(theFactorizationScaleFreeze,GeV) << theShowerApproximation << theDoCompensate << theFreezeGrid << theDetuning << theEventReweight << theSplittingReweight << ounit(maxPt,GeV) << ounit(muPt,GeV)<< theMergingHelper << theColouredOffShellInShower << theInputColouredOffShellInShower << _rearrange << _dipmax << _diplong << _rearrangeNEmissions << theZBoundaries; } void DipoleShowerHandler::persistentInput(PersistentIStream & is, int) { is >> kernels >> theEvolutionOrdering >> constituentReshuffler >> intrinsicPtGenerator >> theGlobalAlphaS >> chainOrderVetoScales >> nEmissions >> discardNoEmissions >> firstMCatNLOEmission >> thePowhegDecayEmission >> realignmentScheme >> verbosity >> printEvent >> iunit(theRenormalizationScaleFreeze,GeV) >> iunit(theFactorizationScaleFreeze,GeV) >> theShowerApproximation >> theDoCompensate >> theFreezeGrid >> theDetuning >> theEventReweight >> theSplittingReweight >> iunit(maxPt,GeV) >> iunit(muPt,GeV)>>theMergingHelper >> theColouredOffShellInShower >> theInputColouredOffShellInShower >> _rearrange >> _dipmax >> _diplong >> _rearrangeNEmissions >> theZBoundaries; } ClassDescription DipoleShowerHandler::initDipoleShowerHandler; // Definition of the static class description member. void DipoleShowerHandler::Init() { static ClassDocumentation documentation ("The DipoleShowerHandler class manages the showering using " "the dipole shower algorithm.", "The shower evolution was performed using the algorithm described in " "\\cite{Platzer:2009jq} and \\cite{Platzer:2011bc}.", "%\\cite{Platzer:2009jq}\n" "\\bibitem{Platzer:2009jq}\n" "S.~Platzer and S.~Gieseke,\n" "``Coherent Parton Showers with Local Recoils,''\n" " JHEP {\\bf 1101}, 024 (2011)\n" "arXiv:0909.5593 [hep-ph].\n" "%%CITATION = ARXIV:0909.5593;%%\n" "%\\cite{Platzer:2011bc}\n" "\\bibitem{Platzer:2011bc}\n" "S.~Platzer and S.~Gieseke,\n" "``Dipole Showers and Automated NLO Matching in Herwig,''\n" "arXiv:1109.6256 [hep-ph].\n" "%%CITATION = ARXIV:1109.6256;%%"); static RefVector interfaceKernels ("Kernels", "Set the splitting kernels to be used by the dipole shower.", &DipoleShowerHandler::kernels, -1, false, false, true, false, false); static Reference interfaceEvolutionOrdering ("EvolutionOrdering", "Set the evolution ordering to be used.", &DipoleShowerHandler::theEvolutionOrdering, false, false, true, false, false); static Reference interfaceConstituentReshuffler ("ConstituentReshuffler", "The object to be used to reshuffle partons to their constitutent mass shells.", &DipoleShowerHandler::constituentReshuffler, false, false, true, true, false); static Reference interfaceIntrinsicPtGenerator ("IntrinsicPtGenerator", "Set the object in charge to generate intrinsic pt for incoming partons.", &DipoleShowerHandler::intrinsicPtGenerator, false, false, true, true, false); static Reference interfaceGlobalAlphaS ("GlobalAlphaS", "Set a global strong coupling for all splitting kernels.", &DipoleShowerHandler::theGlobalAlphaS, false, false, true, true, false); static Switch interfaceRealignmentScheme ("RealignmentScheme", "The realignment scheme to use.", &DipoleShowerHandler::realignmentScheme, 0, false, false); static SwitchOption interfaceRealignmentSchemePreserveRapidity (interfaceRealignmentScheme, "PreserveRapidity", "Preserve the rapidity of non-coloured outgoing system.", 0); static SwitchOption interfaceRealignmentSchemeEvolutionFractions (interfaceRealignmentScheme, "EvolutionFractions", "Use momentum fractions as generated by the evolution.", 1); static SwitchOption interfaceRealignmentSchemeCollisionFrame (interfaceRealignmentScheme, "CollisionFrame", "Determine realignment from collision frame.", 2); static Switch interfaceChainOrderVetoScales ("ChainOrderVetoScales", "[experimental] Switch the chain ordering for veto scales on or off.", &DipoleShowerHandler::chainOrderVetoScales, true, false, false); static SwitchOption interfaceChainOrderVetoScalesYes (interfaceChainOrderVetoScales, "Yes", "Switch on chain ordering for veto scales.", true); static SwitchOption interfaceChainOrderVetoScalesNo (interfaceChainOrderVetoScales, "No", "Switch off chain ordering for veto scales.", false); interfaceChainOrderVetoScales.rank(-1); static Parameter interfaceNEmissions ("NEmissions", "[debug option] Limit the number of emissions to be generated. Zero does not limit the number of emissions.", &DipoleShowerHandler::nEmissions, 0, 0, 0, false, false, Interface::lowerlim); interfaceNEmissions.rank(-1); static Switch interfaceDiscardNoEmissions ("DiscardNoEmissions", "[debug option] Discard events without radiation.", &DipoleShowerHandler::discardNoEmissions, false, false, false); static SwitchOption interfaceDiscardNoEmissionsYes (interfaceDiscardNoEmissions, "Yes", "Discard events without radiation.", true); static SwitchOption interfaceDiscardNoEmissionsNo (interfaceDiscardNoEmissions, "No", "Do not discard events without radiation.", false); interfaceDiscardNoEmissions.rank(-1); static Switch interfaceFirstMCatNLOEmission ("FirstMCatNLOEmission", "[debug option] Only perform the first MC@NLO emission.", &DipoleShowerHandler::firstMCatNLOEmission, false, false, false); static SwitchOption interfaceFirstMCatNLOEmissionYes (interfaceFirstMCatNLOEmission, "Yes", "Perform only the first MC@NLO emission.", true); static SwitchOption interfaceFirstMCatNLOEmissionNo (interfaceFirstMCatNLOEmission, "No", "Produce all emissions.", false); interfaceFirstMCatNLOEmission.rank(-1); static Parameter interfaceVerbosity ("Verbosity", "[debug option] Set the level of debug information provided.", &DipoleShowerHandler::verbosity, 0, 0, 0, false, false, Interface::lowerlim); interfaceVerbosity.rank(-1); static Parameter interfacePrintEvent ("PrintEvent", "[debug option] The number of events for which debugging information should be provided.", &DipoleShowerHandler::printEvent, 0, 0, 0, false, false, Interface::lowerlim); interfacePrintEvent.rank(-1); static Parameter interfaceRenormalizationScaleFreeze ("RenormalizationScaleFreeze", "The freezing scale for the renormalization scale.", &DipoleShowerHandler::theRenormalizationScaleFreeze, GeV, 1.0*GeV, 0.0*GeV, 0*GeV, false, false, Interface::lowerlim); static Parameter interfaceFactorizationScaleFreeze ("FactorizationScaleFreeze", "The freezing scale for the factorization scale.", &DipoleShowerHandler::theFactorizationScaleFreeze, GeV, 2.0*GeV, 0.0*GeV, 0*GeV, false, false, Interface::lowerlim); static Switch interfaceDoCompensate ("DoCompensate", "", &DipoleShowerHandler::theDoCompensate, false, false, false); static SwitchOption interfaceDoCompensateYes (interfaceDoCompensate, "Yes", "", true); static SwitchOption interfaceDoCompensateNo (interfaceDoCompensate, "No", "", false); static Parameter interfaceFreezeGrid ("FreezeGrid", "", &DipoleShowerHandler::theFreezeGrid, 500000, 1, 0, false, false, Interface::lowerlim); static Parameter interfaceDetuning ("Detuning", "A value to detune the overestimate kernel.", &DipoleShowerHandler::theDetuning, 1.0, 1.0, 0, false, false, Interface::lowerlim); static Reference interfaceEventReweight ("EventReweight", "", &DipoleShowerHandler::theEventReweight, false, false, true, true, false); static Reference interfaceSplittingReweight ("SplittingReweight", "Set the splitting reweight.", &DipoleShowerHandler::theSplittingReweight, false, false, true, true, false); static Switch interfacePowhegDecayEmission ("PowhegDecayEmission", "Use Powheg style emission for the decays", &DipoleShowerHandler::thePowhegDecayEmission, true, false, false); static SwitchOption interfacePowhegDecayEmissionYes (interfacePowhegDecayEmission,"Yes","Powheg decay emission on", true); static SwitchOption interfacePowhegDecayEmissionNo (interfacePowhegDecayEmission,"No","Powheg decay emission off", false); static ParVector interfaceOffShellInShower ("OffShellInShower", "PDG codes of the coloured particles that can be off-shell in the process.", &DipoleShowerHandler::theInputColouredOffShellInShower, -1, 0l, -10000000l, 10000000l, false, false, Interface::limited); static Switch interfacerearrange ("Rearrange", "Allow rearranging of dipole chains according to arXiv:1801.06113", &DipoleShowerHandler::_rearrange, false, false, false); static SwitchOption interfacerearrangeYes (interfacerearrange,"Yes","_rearrange on", true); static SwitchOption interfacerearrangeNo (interfacerearrange,"No","_rearrange off", false); static Parameter interfacedipmax ("DipMax", "Allow rearrangment of color chains with ME including dipmax dipoles.", &DipoleShowerHandler::_dipmax, 0, 0, 0, false, false, Interface::lowerlim); static Parameter interfacediplong ("DipLong", "Dipole chains with more than dipmax dipoles are treated as long. \ diplong=3 rearranges these chains with eeuugg MEs, \ diplong=4 rearranges these chains with eeuuggg MEs (slower), \ diplong=5 rearranges these chains with eeuugggg MEs (slow).\ Note: Numerically there is no difference between the options. ", &DipoleShowerHandler::_diplong, 0, 0, 0, false, false, Interface::lowerlim); static Parameter interfacedcorrectNemissions ("RearrangeNEmissions", "Allow rearrangment of color chains up to the nth emission.", &DipoleShowerHandler::_rearrangeNEmissions, 0, 0, 0, false, false, Interface::lowerlim); }