diff --git a/Looptools/util/ffcxs3.F b/Looptools/util/ffcxs3.F --- a/Looptools/util/ffcxs3.F +++ b/Looptools/util/ffcxs3.F @@ -1,654 +1,656 @@ #include "externals.h" #include "types.h" *###[ ffcxs3: subroutine ffcxs3(cs3,ipi12,y,z,dyz,d2yzz,dy2z,xpi,piDpj,ii,ns, + isoort,ier) ***#[*comment:*********************************************************** * * * calculates the s3 as defined in appendix b. * * (ip = ii+3, is1 = ii, is2 = ii+1) * * * * log( xk*y^2 + (-xk+xm1-xm2)*y + xm2 - i*eps ) * * /1 - log( ... ) |y=yi * * s3 = \ dy -------------------------------------------------- * * /0 y - yi * * * * = r(yi,y-,+) + r(yi,y+,-) * * * * with y+- the roots of the argument of the logarithm. * * the sign of the argument to the logarithms in r is passed * * in ieps * * * * input: y(4),z(4) (real) roots in form (z-,z+,1-z-,1-z+) * * dyz(2,2),d2yzz, (real) y() - z(), y+ - z- - z+ * * dy2z(4) (real) y() - 2z() * * xpi (real(ns)) p(i).p(i) (B&D metric) i=1,3 * * m(i)^2 = si.si i=4,6 * * ii (integer) xk = xpi(ii+3) etc * * ns (integer) size of arrays * * isoort (integer) returns kind of action taken * * cs3 (complex)(20) assumed zero. * * ccy (complex)(3) if i0 != 0: complex y * * * * output: cs3 (complex) mod factors pi^2/12, in array * * ipi12 (integer) these factors * * ier (integer) 0=ok 1=inaccurate 2=error * * * * calls: ffcrr,ffcxr,real/dble,ToComplex,log,ffadd1,ffadd2,ffadd3 * * * ***#]*comment:*********************************************************** * #[ declarations: implicit none * * arguments: * integer ipi12(2),ii,ns,isoort(2),ier ComplexType cs3(20) RealType y(4),z(4),dyz(2,2),d2yzz,dy2z(4), + xpi(ns),piDpj(ns,ns) * * local variables: * integer i,ip,ieps(2) RealType yy,yy1,zz,zz1,dyyzz,xdilog,xlog,x00(3) logical ld2yzz * * common blocks * #include "ff.h" * * #] declarations: * #[ get counters: ip = ii+3 if ( isoort(2) .ne. 0 ) then if ( (z(2).gt.z(1) .or. z(1).eq.z(2) .and. z(4).lt.z(3) ) + .eqv. xpi(ip) .gt. 0 ) then ieps(1) = +1 ieps(2) = -1 else ieps(1) = -1 ieps(2) = +1 endif else if ( piDpj(ip,ii) .gt. 0 ) then ieps(1) = +1 else ieps(1) = -1 endif endif * #] get counters: * #[ special case |z| >> |y|: if ( xpi(ip).lt.0 .and. max(abs(y(2)),abs(y(4))) .lt. + xloss*min(abs(z(1)), abs(z(2)))/2 ) then * * we will obtain cancellations of the type Li_2(x) + Li_2(-x) * with x small. * yy = dyz(2,1)/d2yzz yy1 = dyz(2,2)/d2yzz if ( y(2) .eq. 0 ) goto 10 zz = z(2)*yy/y(2) zz1 = 1-zz dyyzz = dyz(2,2)*yy/y(2) call ffcxr(cs3(1),ipi12(1),yy,yy1,zz,zz1,dyyzz,.FALSE., + 0D0,0D0,0D0,.FALSE.,x00,0,ier) 10 continue if ( y(4) .eq. 0 ) goto 30 zz = yy*z(4)/y(4) zz1 = 1-zz dyyzz = -yy*dyz(2,2)/y(4) call ffcxr(cs3(8),ipi12(2),yy,yy1,zz,zz1,dyyzz,.FALSE., + 0D0,0D0,0D0,.FALSE.,x00,0,ier) do 20 i=8,14 - 20 cs3(i) = -cs3(i) + cs3(i) = -cs3(i) + 20 continue 30 continue * And now the remaining Li_2(x^2) terms call ffxli2(xdilog,xlog,(y(2)/dyz(2,1))**2,ier) cs3(15) = +xdilog/2 call ffxli2(xdilog,xlog,(y(4)/dyz(2,1))**2,ier) cs3(16) = -xdilog/2 goto 900 endif * #] special case |z| >> |y|: * #[ normal: if ( xpi(ip) .eq. 0 ) then ld2yzz = .FALSE. else ld2yzz = .TRUE. endif if ( isoort(1) .ne. 0 ) call ffcxr(cs3(1),ipi12(1),y(2),y(4), + z(1),z(3),dyz(2,1),ld2yzz,d2yzz,z(2),z(4),.TRUE.,dy2z(1), + ieps(1),ier) if ( isoort(2) .ne. 0 ) then if ( mod(isoort(2),10) .eq. 2 ) then * both roots are equal: multiply by 2 do 60 i=1,7 cs3(i) = 2*Re(cs3(i)) 60 continue ipi12(1) = 2*ipi12(1) else call ffcxr(cs3(8),ipi12(2),y(2),y(4),z(2),z(4),dyz(2,2), + ld2yzz,d2yzz,z(1),z(3),.TRUE.,dy2z(2),ieps(2),ier) endif endif * * #] normal: 900 continue *###] ffcxs3: end *###[ ffcs3: subroutine ffcs3(cs3,ipi12,cy,cz,cdyz,cd2yzz,cpi,cpiDpj,ii,ns, + isoort,ier) ***#[*comment:*********************************************************** * * * calculates the s3 as defined in appendix b. * * * * log( cpi(ii+3)*y^2 + (cpi(ii+3)+cpi(ii)-cpi(ii+1))*y * * /1 + cpi(ii+1)) - log( ... ) |y=cyi * * s3 = \ dy ---------------------------------------------------- * * /0 y - cyi * * * * = r(cyi,cy+) + r(cyi,cy-) + ( eta(-cy-,-cy+) - * * eta(1-cy-,1-cy+) - eta(...) )*log(1-1/cyi) * * * * with y+- the roots of the argument of the logarithm. * * * * input: cy(4) (complex) cy(1)=y^-,cy(2)=y^+,cy(i+2)=1-cy(1) * * cz(4) (complex) cz(1)=z^-,cz(2)=z^+,cz(i+2)=1-cz(1) * * cpi(6) (complex) masses & momenta (B&D) * * ii (integer) position of cp,cma,cmb in cpi * * ns (integer) size of arrays * * isoort(2)(integer) returns the kind of action taken * * cs3 (complex)(14) assumed zero. * * * * output: cs3 (complex) mod factors ipi12 * * ipi12(2) (integer) these factors * * ier (integer) 0=ok, 1=numerical problems, 2=error * * * * calls: ffcrr,Im,Re,zfflog * * * ***#]*comment:*********************************************************** * #[ declarations: implicit none * * arguments: * integer ipi12(2),ii,ns,isoort(2),ier ComplexType cs3(20),cpi(ns),cpiDpj(ns,ns) ComplexType cy(4),cz(4),cdyz(2,2),cd2yzz * * local variables: * integer i,ip,ieps(2),ieps0,ni(4),ntot logical ld2yzz ComplexType c,zdilog,zlog,cyy,cyy1,czz,czz1,cdyyzz RealType absc,y,y1,z,z1,dyz,d2yzz,zz,zz1, + x00(3),sprec * * common blocks * #include "ff.h" * * statement function * absc(c) = abs(Re(c)) + abs(Im(c)) * #] declarations: * #[ get ieps: ip = ii+3 call ffieps(ieps,cz(1),cpi(ip),cpiDpj(ip,ii),isoort) * #] get ieps: * #[ special case |cz| >> |cy|: if ( isoort(2) .ne. 0 .and. max(absc(cy(2)),absc(cy(4))) .lt. + xloss*min(absc(cz(1)),absc(cz(2)))/2 ) then * * we will obtain cancellations of the type Li_2(x) + Li_2(-x) * with x small. * cyy = cdyz(2,1)/cd2yzz cyy1 = cdyz(2,2)/cd2yzz if ( absc(cy(2)) .lt. xclogm ) then if ( Im(cy(2)) .eq. 0 .and. abs(Re(cy(2))) .gt. + xalogm ) then czz = cz(2)*cyy*ToComplex(1/Re(cy(2))) cdyyzz = cyy*cdyz(2,2)*ToComplex(1/Re(cy(2))) elseif ( cy(2) .eq. 0 .and. cz(2) .ne. 0 .and. cyy + .ne. 0 ) then * the answer IS zero goto 30 endif else czz = cz(2)*cyy/cy(2) cdyyzz = cyy*cdyz(2,2)/cy(2) endif czz1 = 1-czz if ( isoort(1) .eq. -10 ) then * no eta terms. ieps0 = 99 else * do not know the im part ieps0 = 0 endif call ffcrr(cs3(1),ipi12(1),cyy,cyy1,czz,czz1,cdyyzz,.FALSE., + czero,czero,czero,-1,ieps0,ier) 30 continue if ( absc(cy(4)) .lt. xclogm ) then if ( Im(cy(4)) .eq. 0 .and. abs(Re(cy(4))) .gt. + xalogm ) then czz = cz(4)*cyy*ToComplex(1/Re(cy(4))) cdyyzz = -cyy*cdyz(2,2)*ToComplex(1/Re(cy(4))) elseif ( cy(4) .eq. 0 .and. cz(4) .ne. 0 .and. cyy + .ne. 0 ) then * the answer IS zero goto 50 endif else czz = cz(4)*cyy/cy(4) cdyyzz = -cyy*cdyz(2,2)/cy(4) endif czz1 = 1-czz call ffcrr(cs3(8),ipi12(2),cyy,cyy1,czz,czz1,cdyyzz,.FALSE., + czero,czero,czero,-1,ieps0,ier) do 40 i=8,14 cs3(i) = -cs3(i) 40 continue 50 continue * * And now the remaining Li_2(x^2) terms * stupid Gould NP1 * c = cy(2)*cy(2)/(cdyz(2,1)*cdyz(2,1)) call ffzli2(zdilog,zlog,c,ier) cs3(15) = +zdilog/2 * stupid Gould NP1 c = cy(4)*cy(4)/(cdyz(2,1)*cdyz(2,1)) call ffzli2(zdilog,zlog,c,ier) cs3(16) = -zdilog/2 goto 900 endif * #] special case |cz| >> |cy|: * #[ normal: if ( isoort(2) .eq. 0 ) then ld2yzz = .FALSE. else ld2yzz = .TRUE. endif if ( isoort(1) .eq. 0 ) then * do nothing elseif ( mod(isoort(1),10).eq.0 .or. mod(isoort(1),10).eq.-1 + .or. mod(isoort(1),10).eq.-3 ) then call ffcrr(cs3(1),ipi12(1),cy(2),cy(4),cz(1),cz(3), + cdyz(2,1),ld2yzz,cd2yzz,cz(2),cz(4),isoort(1), + ieps(1),ier) elseif ( mod(isoort(1),10) .eq. -5 .or. mod(isoort(1),10) .eq. + -6 ) then y = Re(cy(2)) y1 = Re(cy(4)) z = Re(cz(1)) z1 = Re(cz(3)) dyz = Re(cdyz(2,1)) d2yzz = Re(cd2yzz) zz = Re(cz(2)) zz1 = Re(cz(4)) sprec = precx precx = precc call ffcxr(cs3(1),ipi12(1),y,y1,z,z1,dyz,ld2yzz,d2yzz,zz,zz1 + ,.FALSE.,x00,ieps(1),ier) precx = sprec else call fferr(12,ier) endif if ( isoort(2) .eq. 0 ) then * do nothing elseif ( mod(isoort(2),5) .eq. 0 ) then - do 100 i=1,7 - 100 cs3(i) = 2*Re(cs3(i)) + do 100 i=1,7 + cs3(i) = 2*Re(cs3(i)) + 100 continue ipi12(1) = 2*ipi12(1) elseif ( mod(isoort(2),10).eq.-1 .or. mod(isoort(1),10).eq.-3 ) + then call ffcrr(cs3(8),ipi12(2),cy(2),cy(4),cz(2),cz(4), + cdyz(2,2),ld2yzz,cd2yzz,cz(1),cz(3),isoort(2), + ieps(2),ier) elseif ( mod(isoort(2),10) .eq. -6 ) then y = Re(cy(2)) y1 = Re(cy(4)) z = Re(cz(2)) z1 = Re(cz(4)) dyz = Re(cdyz(2,2)) d2yzz = Re(cd2yzz) zz = Re(cz(1)) zz1 = Re(cz(3)) sprec = precx precx = precc call ffcxr(cs3(8),ipi12(2),y,y1,z,z1,dyz,ld2yzz,d2yzz,zz,zz1 + ,.FALSE.,x00,ieps(2),ier) precx = sprec else call fferr(13,ier) endif * #] normal: * #[ eta's: if ( mod(isoort(1),10).eq.-5 .or. mod(isoort(1),10).eq.-6 ) + then if ( mod(isoort(2),10).ne.-5 .and. mod(isoort(1),10).ne.-6 + ) then print *,'ffcxs3: error: I assumed both would be real!' ier = ier + 50 endif * we called ffcxr - no eta's elseif ( Im(cpi(ip)).eq.0 ) then call ffgeta(ni,cz(1),cdyz(1,1), + cpi(ip),cpiDpj(ii,ip),ieps,isoort,ier) ntot = ni(1) + ni(2) + ni(3) + ni(4) if ( ntot .ne. 0 ) call ffclgy(cs3(15),ipi12(2),ntot, + cy(1),cz(1),cd2yzz,ier) else * * cpi(ip) is really complex (occurs in transformed * 4pointfunction) * print *,'THIS PART IS NOT READY ', + 'and should not be reached' c stop endif * #] eta's: 900 continue *###] ffcs3: end *###[ ffclgy: subroutine ffclgy(cs3,ipi12,ntot,cy,cz,cd2yzz,ier) ***#[*comment:*********************************************************** * * * calculates the the difference of two S's with cy(3,4),cz(3,4), * * cy(4)cz(3)-cy(3)cz(4) given. Note the difference with ffdcs4, * * in which the cy's are the same and only the cz's different. * * Here both can be different. Also we skip an intermediat * * level. * * * * input: cy(4) (complex) cy,1-cy in S with s3,s4 * * cz(4) (complex) cz,1-cz in S with s3,s4 * * cdyz(2,2) (complex) cy - cz * * cd2yzz (complex) 2*cy - cz+ - cz- * * cdyzzy(4) (complex) cy(i,4)*cz(i,4)-cy(i,3)*cz(i,4) * * cpiDpj(6,6) (complex) usual * * cs3 (complex) assumed zero. * * * * output: cs3 (complex) mod factors pi^2/12, in array * * ipi12 (integer) these factors * * isoort (integer) returns kind of action taken * * ier (integer) number of digits lost * * * ***#]*comment:*********************************************************** * #[ declarations: implicit none * * arguments * ComplexType cs3 ComplexType cy(4),cz(4),cd2yzz integer ipi12,ntot,ier * * local variables * integer ipi ComplexType c,cc,clogy,c2y1,zfflog,zfflo1,csum RealType absc external zfflog,zfflo1 * * common blocks * #include "ff.h" * * statement function * absc(c) = abs(Re(c)) + abs(Im(c)) * #] declarations: * #[ calculations: ipi = 0 if ( 1 .lt. xloss*absc(cy(2)) ) then clogy = zfflo1(1/cy(2),ier) else if ( absc(cy(2)) .lt. xclogm .or. absc(cy(4)) .lt. xclogm ) + then if ( ntot .ne. 0 ) call fferr(15,ier) clogy = 0 else c = -cy(4)/cy(2) if ( Re(c) .gt. -abs(Im(c)) ) then clogy = zfflog(c,0,czero,ier) else * take out the factor 2*pi^2 cc = c+1 if ( absc(cc) .lt. xloss ) then c2y1 = -cd2yzz - cz(1) + cz(4) if ( absc(c2y1) .lt. xloss*max(absc(cz(1)), + absc(cz(4))) ) then c2y1 = -cd2yzz - cz(2) + cz(3) endif csum = -c2y1/cy(2) clogy = zfflo1(csum,ier) else csum = 0 clogy = zfflog(-c,0,czero,ier) endif if ( Im(c) .lt. -precc*absc(c) .or. + Im(csum) .lt. -precc*absc(csum) ) then ipi = -1 elseif ( Im(c) .gt. precc*absc(c) .or. + Im(csum) .gt. precc*absc(csum) ) then ipi = +1 else call fferr(51,ier) ipi = 0 endif endif endif endif cs3 = cs3 + ntot*c2ipi*clogy if ( ipi .ne. 0 ) then ipi12 = ipi12 - 24*ntot*ipi endif * #] calculations: *###] ffclgy: end *###[ ffieps: subroutine ffieps(ieps,cz,cp,cpDs,isoort) ***#[*comment:*********************************************************** * * * Get the ieps prescription in such a way that it is compatible * * with the imaginary part of cz if non-zero, compatible with the * * real case if zero. * * * * Input: cz complex(4) the roots z-,z+,1-z-,1-z+ * * cp complex p^2 * * cpDs complex p.s * * isoort integer(2) which type of Ri * * * * Output: ieps integer(2) z -> z-ieps*i*epsilon * * will give correct im part * * * ***#]*comment:*********************************************************** * #[ declarations: implicit none * * arguments: * integer ieps(2),isoort(2) ComplexType cp,cpDs,cz(4) * * #] declarations: * #[ work: if ( Im(cp) .ne. 0 ) then * do not calculate ANY eta terms, we'll do them ourselves. ieps(1) = 99 ieps(2) = 99 elseif ( isoort(2) .ne. 0 ) then if ( Im(cz(1)) .lt. 0 ) then ieps(1) = +1 if ( Im(cz(2)) .lt. 0 ) then ieps(2) = +1 else ieps(2) = -1 endif elseif ( Im(cz(1)) .gt. 0 ) then ieps(1) = -1 if ( Im(cz(2)) .le. 0 ) then ieps(2) = +1 else ieps(2) = -1 endif else if ( Im(cz(2)) .lt. 0 ) then ieps(1) = -1 ieps(2) = +1 elseif ( Im(cz(2)) .gt. 0 ) then ieps(1) = +1 ieps(2) = -1 else if ( (Re(cz(2)).gt.Re(cz(1)) + .or. (Re(cz(1)).eq.Re(cz(2)) + .and. Re(cz(4)).lt.Re(cz(3))) + ) .eqv. Re(cp).gt.0 ) then ieps(1) = +1 ieps(2) = -1 else ieps(1) = -1 ieps(2) = +1 endif endif endif else if ( Im(cz(1)) .lt. 0 ) then ieps(1) = +1 elseif ( Im(cz(1)) .gt. 0 ) then ieps(1) = -1 elseif ( Re(cpDs) .gt. 0 ) then ieps(1) = +1 else ieps(1) = -1 endif ieps(2) = -9999 endif * #] work: *###] ffieps: end *###[ ffgeta: subroutine ffgeta(ni,cz,cdyz,cp,cpDs,ieps,isoort,ier) ***#[*comment:*********************************************************** * * * Get the eta terms which arise from splitting up * * log(p2(x-z-)(x-z+)) - log(p2(y-z-)(y-z+)) * * * * Input: cz complex(4) the roots z-,z+,1-z-,1-z+ * * cdyz complex(2,2) y-z * * cd2yzz complex(2) 2y-(z-)-(z+) * * cp complex p^2 * * cpDs complex p.s * * ieps integer(2) the assumed im part if Im(z)=0 * * isoort integer(2) which type of Ri * * * * Output: ni integer(4) eta()/(2*pi*i) * * * ***#]*comment:*********************************************************** * #[ declarations: implicit none * * arguments: * integer ni(4),ieps(2),isoort(2),ier ComplexType cp,cpDs,cz(4),cdyz(2,2) * * local variables * integer i,nffeta,nffet1 ComplexType cmip external nffeta,nffet1 * * common * #include "ff.h" * * #] declarations: * #[ complex masses or imaginary roots: * * only complex because of complex roots in y or z * [checked and in agreement with ieps definition 23-sep-1991] * * isoort = +1: y is real, z is real * isoort = -1-n*10: y is complex, possibly z as well * isoort = -3-n*10: y,z complex, (y-z-)*(y-z+) real * isoort = 0: y is complex, one z root only * isoort = -10-n*10: y is real, z is complex * isoort = -5,6-n*10: y,z real * if ( isoort(1) .gt. 0 ) then * * really a real case * ni(1) = 0 ni(2) = 0 ni(3) = 0 ni(4) = 0 elseif ( mod(isoort(1),10) .ne. 0 .and. isoort(2) .ne. 0 ) then cmip = ToComplex(0D0,-Re(cp)) * * ni(1) = eta(p2,(x-z-)(x-z+)) = 0 by definition (see ni(3)) * ni(2) = eta(x-z-,x-z+) * ni(1) = 0 if ( ieps(1) .gt. 0 .neqv. ieps(2) .gt. 0 ) then ni(2) = 0 else ni(2) = nffet1(-cz(1),-cz(2),cmip,ier) if ( cz(3).ne.0 .and. cz(4).ne.0 ) then i = nffet1(cz(3),cz(4),cmip,ier) if ( i .ne. ni(2) ) call fferr(53,ier) endif endif * * ni(3) compensates for whatever convention we chose in ni(1) * ni(4) = -eta(y-z-,y-z+) * if ( mod(isoort(1),10).eq.-3 ) then * follow the i*epsilon prescription as (y-z-)(y-z+) real ni(3) = 0 ni(4) = -nffet1(cdyz(2,1),cdyz(2,2),cmip,ier) else if ( Re(cp) .lt. 0 .and. Im(cdyz(2,1)* + cdyz(2,2)) .lt. 0 ) then ni(3) = -1 else ni(3) = 0 endif ni(4) = -nffeta(cdyz(2,1),cdyz(2,2),ier) endif elseif ( (mod(isoort(1),10).eq.-1 .or. mod(isoort(1),10).eq.-3) + .and. isoort(2) .eq. 0 ) then ni(1) = 0 if ( Im(cz(1)) .ne. 0 ) then ni(2) = nffet1(-cpDs,-cz(1),ToComplex(Re(0), + Re(-1)),ier) else ni(2) = nffet1(-cpDs,ToComplex(Re(0),Re(1)), + ToComplex(Re(0),Re(-1)),ier) endif ni(3) = 0 ni(4) = -nffeta(-cpDs,cdyz(2,1),ier) else ni(1) = 0 ni(2) = 0 ni(3) = 0 ni(4) = 0 endif * #] complex masses or imaginary roots: *###] ffgeta: end diff --git a/Looptools/util/ffcxs4.F b/Looptools/util/ffcxs4.F --- a/Looptools/util/ffcxs4.F +++ b/Looptools/util/ffcxs4.F @@ -1,778 +1,784 @@ #include "externals.h" #include "types.h" * $Id: ffcxs4.f,v 1.3 1995/10/17 06:55:09 gj Exp $ * $Log: ffcxs4.f,v $ c Revision 1.3 1995/10/17 06:55:09 gj c Fixed ieps error in ffdcrr (ffcxs4.f), added real case in ffcrr, debugging c info in ffxd0, and warned against remaining errors for del2=0 in ffrot4 c (ffxd0h.f) c c Revision 1.2 1995/10/06 09:17:22 gj c Found stupid typo in ffxc0p which caused the result to be off by pi^2/3 in c some equal-mass cases. Added checks to ffcxs4.f ffcrr.f. c *###[ ffcxs4: subroutine ffcxs4(cs3,ipi12,w,y,z,dwy,dwz,dyz,d2yww,d2yzz, + xpi,piDpj,ii,ns,isoort,ier) ***#[*comment:*********************************************************** * * * Calculate the 8 Spence functions = 4 R's = 2 dR's * * * * * * * ***#]*comment:*********************************************************** * #[ declarations: implicit none * * arguments: * integer ipi12(4),ii,ns,isoort(4),ier ComplexType cs3(40) RealType w(4),y(4),z(4),dwy(2,2),dwz(2,2),dyz(2,2), + d2yww,d2yzz,xpi(ns),piDpj(ns,ns),x00(3) * * local variables: * integer iepz(2),iepw(2) logical ld2yzz,ld2yww * * common blocks * #include "ff.h" * #] declarations: * #[ groundwork: if ( isoort(2) .eq. 0 ) then ld2yzz = .FALSE. else ld2yzz = .TRUE. endif if ( isoort(4) .eq. 0 ) then ld2yww = .FALSE. else ld2yww = .TRUE. endif if ( isoort(2) .ne. 0 ) then if ( z(2) .gt. z(1) .eqv. xpi(ii+3) .gt. 0 ) then iepz(1) = +1 iepz(2) = -1 else iepz(1) = -1 iepz(2) = +1 endif else print *,'ffcxs4: error: untested algorithm' if ( piDpj(ii,ii+3) .gt. 0 ) then iepz(1) = +1 else iepz(1) = -1 endif endif if ( isoort(4) .ne. 0 ) then if ( w(2) .gt. w(1) .eqv. xpi(5) .gt. 0 ) then iepw(1) = 1 iepw(2) = -1 else iepw(1) = -1 iepw(2) = 1 endif else print *,'ffcxs4: error: untested algorithm' if ( piDpj(2,5) .gt. 0 ) then iepw(1) = +1 else iepw(1) = -1 endif endif * #] groundwork: * #[ zm and wp: if ( isoort(4) .eq. 0 ) then call ffcxr(cs3(1),ipi12(1),y(2),y(4),z(1),z(3),dyz(2,1), + ld2yzz,d2yzz,z(2),z(4),.FALSE.,x00,iepz(1),ier) else if ( .not. ( dwz(2,1).eq.0 .and. iepz(1).eq.iepw(2) ) ) + call ffdcxr(cs3( 1),ipi12(1),y(2),y(4),z(1),z(3), + z(2),z(4),d2yzz,w(2),w(4),w(1),w(3),d2yww, + dyz(2,1),dwy(2,2),dwz(2,1),iepz(1),iepw(2),ier) endif * #] zm and wp: * #[ zp and wm: if ( isoort(2) .eq. 0 ) then call ffcxr(cs3(1),ipi12(1),y(2),y(4),w(1),w(3),-dwy(1,2), + ld2yww,d2yww,w(2),w(4),.FALSE.,x00,iepw(1),ier) else if ( .not. ( dwz(1,2).eq.0 .and. iepz(2).eq.iepw(1) ) ) + call ffdcxr(cs3(21),ipi12(3),y(2),y(4),z(2),z(4), + z(1),z(3),d2yzz,w(1),w(3),w(2),w(4),d2yww, + dyz(2,2),dwy(1,2),dwz(1,2),iepz(2),iepw(1),ier) endif * #] zp and wm: *###] ffcxs4: end *###[ ffcs4: subroutine ffcs4(cs3,ipi12,cw,cy,cz,cdwy,cdwz,cdyz,cd2yww,cd2yzz + ,cpi,cpiDpj,cp2p,ii,ns,isoort,ier) ***#[*comment:*********************************************************** * * * Calculate the 8 Spence functions = 4 R's = 2 dR's * * * * * * * ***#]*comment:*********************************************************** * #[ declarations: implicit none * * arguments: * integer ipi12(4),ii,ns,isoort(4),ier ComplexType cs3(40) ComplexType cw(4),cy(4),cz(4),cdwy(2,2),cdwz(2,2),cdyz(2,2) ComplexType cd2yww,cd2yzz,cpi(ns),cp2p,cpiDpj(ns,ns) * * local variables: * logical ld2yzz,ld2yww integer i,j,ip,iepz(2),iepw(2),nz(4),nw(4),ntot,i2pi ComplexType c,cc,clogy,c2y1,cdyw(2,2) ComplexType zfflo1,zfflog RealType absc external zfflo1,zfflog * * common blocks * #include "ff.h" * * statement function * absc(c) = abs(Re(c)) + abs(Im(c)) * #] declarations: * #[ get counters: ip = ii+3 if ( isoort(2) .eq. 0 ) then ld2yzz = .FALSE. else ld2yzz = .TRUE. endif if ( isoort(4) .eq. 0 ) then ld2yww = .FALSE. else ld2yww = .TRUE. endif call ffieps(iepz,cz,cpi(ip),cpiDpj(ip,ii),isoort) call ffieps(iepw,cw,cp2p,cpiDpj(ip,ii),isoort(3)) if ( isoort(4) .eq. 0 ) then print *,'ffcs4: error: case not implemented' ier = ier + 50 endif * #] get counters: * #[ R's: if ( isoort(4) .eq. 0 ) then call ffcrr(cs3(1),ipi12(1),cy(2),cy(4),cz(1),cz(3),cdyz(2,1) + ,ld2yzz,cd2yzz,cz(2),cz(4),isoort(4),iepz(1),ier) else if ( .not. ( cdwz(2,1).eq.0 .and. iepz(1).eq.iepw(2) ) ) + call ffdcrr(cs3( 1),ipi12(1),cy(2),cy(4),cz(1),cz(3),cz(2), + cz(4),cd2yzz,cw(2),cw(4),cw(1),cw(3),cd2yww,cdyz(2,1), + cdwy(2,2),cdwz(2,1),isoort(4),iepz(1),iepw(2),ier) endif if ( isoort(2) .eq. 0 ) then call ffcrr(cs3(1),ipi12(1),cy(2),cy(4),cw(1),cw(3),-cdwy(1,2 + ),ld2yww,cd2yww,cw(2),cw(4),isoort(2),iepw(1),ier) else if ( .not. ( cdwz(1,2).eq.0 .and. iepz(2).eq.iepw(1) ) ) + call ffdcrr(cs3(21),ipi12(3),cy(2),cy(4),cz(2),cz(4),cz(1), + cz(3),cd2yzz,cw(1),cw(3),cw(2),cw(4),cd2yww,cdyz(2,2), + cdwy(1,2),cdwz(1,2),iepz(2),isoort(2),iepw(1),ier) endif * #] R's: * #[ eta's: if ( Im(cpi(ip)) .eq. 0 ) then call ffgeta(nz,cz,cdyz, + cpi(ip),cpiDpj(ii,ip),iepz,isoort,ier) do 120 i=1,2 do 110 j=1,2 cdyw(i,j) = cdwy(j,i) 110 continue 120 continue call ffgeta(nw,cw,cdyw, + cp2p,cpiDpj(ii,ip),iepw,isoort(3),ier) else print *,'ffcs4: error: not ready for complex D0 yet' endif ntot = nz(1)+nz(2)+nz(3)+nz(4)-nw(1)-nw(2)-nw(3)-nw(4) if ( ntot .ne. 0 ) then i2pi = 0 if ( 1/absc(cy(2)) .lt. xloss ) then clogy = zfflo1(1/cy(2),ier) else c = -cy(4)/cy(2) if ( Re(c) .gt. -abs(Im(c)) ) then clogy = zfflog(c,0,czero,ier) else * take out the factor 2*pi^2 cc = c+1 if ( absc(cc) .lt. xloss ) then c2y1 = -cd2yzz - cz(1) + cz(4) if ( absc(c2y1) .lt. xloss*max(absc(cz(1)), + absc(cz(4))) ) then c2y1 = -cd2yzz - cz(2) + cz(3) endif clogy = zfflo1(-c2y1/cy(2),ier) else clogy = zfflog(-c,0,czero,ier) endif if ( Im(c) .lt. 0 ) then i2pi = -1 elseif ( Im(c) .gt. 0 ) then i2pi = +1 else call fferr(51,ier) i2pi = 0 endif ipi12(2) = ipi12(2) - ntot*24*i2pi endif endif if ( cs3(40) .ne. 0 ) print *,'ffcs4: error: cs3(40) != 0' cs3(40) = ntot*c2ipi*clogy endif * #] eta's: *###] ffcs4: end *###[ ffdcxr: subroutine ffdcxr(cs3,ipi12,y,y1,z,z1,zp,zp1,d2yzz, + w,w1,wp,wp1,d2yww,dyz,dwy,dwz,iepsz,iepsw,ier) ***#[*comment:*********************************************************** * * * Calculate * * * * R(y,z,iepsz) - R(y,w,iepsw) * * * * Input: * * a = [yzw] (real) see definition * * a1 = 1 - a (real) * * dab = a - b (real) * * ieps[zw] (integer) sign of imaginary part * * of argument logarithm * * cs3(20) (complex) assumed zero * * * * Output: * * cs3(20) (complex) the results, not added * * ipi12(2) (integer) factors pi^2/12 * * * * Calls: ffcxr * * * ***#]*comment:*********************************************************** * #[ declarations: implicit none * * arguments: * integer ipi12(2),iepsz,iepsw,ier ComplexType cs3(20) RealType y,z,w,y1,z1,w1,dyz,dwy,dwz,zp,zp1,d2yzz,wp,wp1, + d2yww * * local variables: * integer i,ieps logical again RealType yy,yy1,zz,zz1,dyyzz,xx1,xx1n,term,tot,d2,d3, + d21,d31,d2n,d3n,d21n1,d31n1,dw,x00(3) ComplexType chulp RealType dfflo1 external dfflo1 * * common blocks * #include "ff.h" * #] declarations: * #[ groundwork: if ( dwz .eq. 0 .and. iepsz .eq. iepsw ) return if ( dyz .eq. 0 ) then call fferr(75,ier) return endif xx1 = y/dyz dw = dwz/dyz if ( xx1 .le. .5D0 .or. xx1 .gt. 2 ) then d2 = 1/y dw = dw*y/w else d2 = 1/z1 endif again = .FALSE. 123 continue * #] groundwork: * #[ trivial case: if ( dw .eq. 0 ) then * #] trivial case: * #[ normal case: elseif ( abs(dw) .gt. xloss .or. again ) then * nothing's the matter call ffcxr(cs3( 1),ipi12(1),y,y1,z,z1,dyz, + .TRUE.,d2yzz,zp,zp1,.FALSE.,x00,iepsz,ier) call ffcxr(cs3(11),ipi12(2),y,y1,w,w1,-dwy, + .TRUE.,d2yww,wp,wp1,.FALSE.,x00,iepsw,ier) do 10 i=11,20 - 10 cs3(i) = -cs3(i) + cs3(i) = -cs3(i) + 10 continue ipi12(2) = -ipi12(2) * #] normal case: * #[ only cancellations in w, not in y: elseif ( abs(d2) .gt. xloss ) then * there are no cancellations the other way: if ( iepsz .ne. iepsw .and. ( y/dyz .gt. 1 .or.-y/dwy .gt. + 1 ) ) then again = .TRUE. goto 123 endif yy = dwy/dwz zz = yy*z/y yy1 = dyz/dwz zz1 = yy1*w/y dyyzz = yy*dyz/y if ( y .lt. 0 ) then ieps = iepsz else ieps = -iepsz endif call ffcxr(cs3( 1),ipi12(1),yy,yy1,zz,zz1,dyyzz,.FALSE., + 0D0,0D0,0D0,.FALSE.,x00,2*ieps,ier) zz = yy*z1/y1 zz1 = yy1*w1/y1 dyyzz = -yy*dyz/y1 if ( y1 .gt. 0 ) then ieps = iepsz else ieps = -iepsz endif call ffcxr(cs3(11),ipi12(2),yy,yy1,zz,zz1,dyyzz,.FALSE., + 0D0,0D0,0D0,.FALSE.,x00,2*ieps,ier) do 20 i=11,20 cs3(i) = -cs3(i) 20 continue ipi12(2) = -ipi12(2) * #] only cancellations in w, not in y: * #[ Hill identity: elseif ( ( 1 .gt. xloss*abs(y) .or. abs(xx1) .gt. xloss ) + .and. ( 1 .gt. xloss*abs(z) .or. abs(z/dyz) .gt. xloss ) + .and. ( 1 .gt. xloss*abs(y) .or. abs(dyz/y) .gt. xloss ) + ) then * do a Hill identity on the y,y-1 direction yy = -y*w1/dwy yy1 = w*y1/dwy zz = -z*w1/dwz zz1 = w*z1/dwz dyyzz = -w*w1*(dyz/(dwy*dwz)) if ( y*dwz .gt. 0 .eqv. (y+dwz) .gt. 0 ) then ieps = 2*iepsw else ieps = -2*iepsw endif call ffcxr(cs3( 1),ipi12(1),yy,yy1,zz,zz1,dyyzz,.FALSE., + 0D0,0D0,0D0,.FALSE.,x00,ieps,ier) yy = w1 yy1 = w zz = -w1*z/dwz zz1 = w*z1/dwz dyyzz = w*w1/dwz call ffcxr(cs3( 9),ipi12(2),yy,yy1,zz,zz1,dyyzz,.FALSE., + 0D0,0D0,0D0,.FALSE.,x00,ieps,ier) do 30 i=9,16 - 30 cs3(i) = -cs3(i) + cs3(i) = -cs3(i) + 30 continue ipi12(2) = -ipi12(2) * the extra logarithms ... if ( 1 .lt. xloss*abs(w) ) then chulp = dfflo1(1/w,ier) elseif ( w1 .lt. 0 .or. w .lt. 0 ) then chulp = log(-w1/w) else chulp = ToComplex(Re(log(w1/w)),Re(-iepsw*pi)) endif cs3(20) = -Re(dfflo1(dwz/dwy,ier))*chulp * #] Hill identity: * #[ Taylor expansion: elseif ( (w.lt.0..or.w1.lt.0) .and. (z.lt.0..or.z1.lt.0) ) then * do a Taylor expansion if ( abs(xx1) .lt. xloss ) then d3 = dwz/dwy xx1n = xx1 d2n = d2 d3n = d3 d21 = 1-d2 d21n1 = 1 d31 = 1-d3 d31n1 = 1 tot = xx1*d2*d3 do 50 i=2,20 xx1n = xx1n*xx1 d21n1 = d21n1*d21 d31n1 = d31n1*d31 d2n = d2n + d2*d21n1 d3n = d3n + d3*d31n1 term = xx1n*d2n*d3n*xn2inv(i) tot = tot + term if ( abs(term) .le. precx*abs(tot) ) goto 51 50 continue 51 continue cs3(1) = tot elseif ( abs(z/dyz) .lt. xloss ) then call ffcxr(cs3( 1),ipi12(1),y,y1,z,z1,dyz, + .TRUE.,d2yzz,zp,zp1,.FALSE.,x00,iepsz,ier) call ffcxr(cs3(11),ipi12(2),y,y1,w,w1,-dwy, + .TRUE.,d2yww,wp,wp1,.FALSE.,x00,iepsw,ier) do 110 i=11,20 - 110 cs3(i) = -cs3(i) + cs3(i) = -cs3(i) + 110 continue else call fferr(22,ier) return endif else call ffcxr(cs3( 1),ipi12(1),y,y1,z,z1,dyz,.FALSE., + 0D0,0D0,0D0,.FALSE.,x00,iepsz,ier) call ffcxr(cs3(11),ipi12(2),y,y1,w,w1,-dwy,.FALSE., + 0D0,0D0,0D0,.FALSE.,x00,iepsw,ier) do 40 i=11,20 - 40 cs3(i) = -cs3(i) + cs3(i) = -cs3(i) + 40 continue ipi12(2) = -ipi12(2) endif * #] Taylor expansion: *###] ffdcxr: end *###[ ffdcrr: subroutine ffdcrr(cs3,ipi12,cy,cy1,cz,cz1,czp,czp1,cd2yzz,cw,cw1 + ,cwp,cwp1,cd2yww,cdyz,cdwy,cdwz,isoort,iepsz,iepsw,ier) ***#[*comment:*********************************************************** * * * Calculate * * * * R(cy,cz,iepsz) - R(cy,cw,iepsw) * * * * Input: * * a = [yzw] (real) see definition * * a1 = 1 - a (real) * * dab = a - b (real) * * ieps[zw] (integer) sign of imaginary part * * of argument logarithm * * cs3(20) (complex) assumed zero * * * * Output: * * cs3(20) (complex) the results, not added * * ipi12(2) (integer) factors pi^2/12 * * * * Calls: ffcrr * * * ***#]*comment:*********************************************************** * #[ declarations: implicit none * * arguments: * integer ipi12(2),isoort,iepsz,iepsw,ier ComplexType cs3(20) ComplexType cy,cz,czp,cw,cwp,cy1,cz1,czp1,cw1,cwp1, + cdyz,cdwy,cdwz,cd2yzz,cd2yww * * local variables: * integer i,ieps,ieps1,ieps2, + nffeta,nffet1,n1,n2,n3,n4,n5,n6 logical ld2yyz ComplexType cyy,cyy1,czz,czz1,cdyyzz,chulp,zfflo1,zfflog, + cc1,cdw,cc1n,cterm,ctot,cd2,cd3, + cd21,cd31,cd2n,cd3n,cd21n1,cd31n1, + cc2,cfactz,cfactw,czzp,czzp1,cd2yyz ComplexType c RealType absc external nffeta,nffet1,zfflo1,zfflog * * common blocks * #include "ff.h" * * statement function * absc(c) = abs(Re(c)) + abs(Im(c)) * #] declarations: * #[ groundwork: if ( cdwz .eq. 0 ) then if ( abs(Im(cz)) .gt. precc*abs(Re(cz)) .or. + iepsz .eq. iepsw ) return if ( Re(cz) .ge. 0 .and. Re(cz1) .ge. 0 ) return call fferr(76,ier) return endif if ( cdyz .eq. 0 ) then call fferr(77,ier) return endif cc1 = cy/cdyz cdw = cdwz/cdyz if ( Re(cc1) .le. .5D0 .or. abs(cc1-1) .gt. 1 ) then cd2 = 1/cy cdw = cdw*cy/cw else cd2 = 1/cz1 endif * #] groundwork: * #[ trivial case: if ( absc(cdw) .eq. 0 ) then * #] trivial case: * #[ normal case: * * if no cancellations are expected OR the imaginary signs differ * and are significant * elseif ( absc(cdw) .gt. xloss .or. (iepsz.ne.iepsw .and. + (Re(cy/cdyz).gt.1 .or. Re(-cy1/cdyz).gt.1) ) ) then * nothing's the matter * special case to avoid bug found 15-oct=1995 if ( iepsz.eq.iepsw ) then if ( Im(cz).eq.0 .and. Im(cz1).eq.0 ) then print *,'ffdcrr: flipping sign iepsz' iepsz = -iepsz elseif ( Im(cw).eq.0 .and. Im(cw1).eq.0 ) then print *,'ffdcrr: flipping sign iepsw' iepsw = -iepsw else print *,'ffdcrr: error: missing eta terms!' ier = ier + 100 endif endif call ffcrr(cs3(1),ipi12(1),cy,cy1,cz,cz1,cdyz,.TRUE., + cd2yzz,czp,czp1,isoort,iepsz,ier) call ffcrr(cs3(8),ipi12(2),cy,cy1,cw,cw1,-cdwy,.TRUE., + cd2yww,cwp,cwp1,isoort,iepsw,ier) do 10 i=8,14 cs3(i) = -cs3(i) 10 continue ipi12(2) = -ipi12(2) * #] normal case: * #[ only cancellations in cw, not in cy: elseif ( absc(cd2) .gt. xloss ) then * there are no cancellations the other way: cyy = cdwy/cdwz czz = cz*cyy/cy cyy1 = cdyz/cdwz czz1 = cyy1*cw/cy cdyyzz = cdyz*cyy/cy if ( Re(cy) .gt. 0 ) then ieps1 = -3*iepsz else ieps1 = +3*iepsz endif * Often 2y-z-z is relevant, but 2*yy-zz-zz is not, solve by * introducing zzp. czzp = czp*cyy/cy cd2yyz = cd2yzz*cyy/cy czzp1 = 1 - czzp if ( absc(czzp1) .lt. xloss ) then * later try more possibilities ld2yyz = .FALSE. else ld2yyz = .TRUE. endif call ffcrr(cs3(1),ipi12(1),cyy,cyy1,czz,czz1,cdyyzz, + ld2yyz,cd2yyz,czzp,czzp1,isoort,ieps1,ier) czz = cyy*cz1/cy1 czz1 = cyy1*cw1/cy1 if ( Re(-cy1) .gt. 0 ) then ieps2 = -3*iepsz else ieps2 = +3*iepsz endif cdyyzz = -cyy*cdyz/cy1 czzp = czp1*cyy/cy1 cd2yyz = -cd2yzz*cyy/cy1 czzp1 = 1 - czzp if ( absc(czzp1) .lt. xloss ) then * later try more possibilities ld2yyz = .FALSE. else ld2yyz = .TRUE. endif call ffcrr(cs3(8),ipi12(2),cyy,cyy1,czz,czz1,cdyyzz, + .TRUE.,cd2yyz,czzp,czzp1,isoort,ieps2,ier) do 20 i=8,14 cs3(i) = -cs3(i) 20 continue ipi12(2) = -ipi12(2) * eta terms (are not calculated in ffcrr as ieps = 3) cfactz = 1/cdyz if ( Im(cz) .eq. 0 ) then if ( Im(cy) .eq. 0 ) then n1 = 0 n2 = 0 else n1 = nffet1(ToComplex(Re(0),Re(iepsz)),cfactz, + -cz*cfactz,ier) n2 = nffet1(ToComplex(Re(0),Re(iepsz)),cfactz, + cz1*cfactz,ier) endif else n1 = nffeta(-cz,cfactz,ier) n2 = nffeta(cz1,cfactz,ier) endif cfactw = -1/cdwy if ( Im(cw) .eq. 0 ) then if ( Im(cy) .eq. 0 ) then n4 = 0 n5 = 0 else n4 = nffet1(ToComplex(Re(0),Re(iepsw)),cfactw, + -cw*cfactw,ier) n5 = nffet1(ToComplex(Re(0),Re(iepsw)),cfactw, + cw1*cfactw,ier) endif else n4 = nffeta(-cw,cfactw,ier) n5 = nffeta(cw1,cfactw,ier) endif * * we assume that cs3(15-17) are not used, this is always true * n3 = 0 n6 = 0 if ( n1.eq.n4 ) then if ( n1.eq.0 ) then * nothing to do else cc1 = cdwz/cdyz if ( absc(cc1) .lt. xloss ) then cs3(15) = n1*c2ipi*zfflo1(cc1,ier) else cc1 = -cdwy/cdyz cs3(15) = n1*c2ipi*zfflog(cc1,0,czero,ier) endif cc1 = cy*cfactz cc2 = cy*cfactw if ( Im(cc1).eq.0 .or. Im(cc2).eq.0 ) then n3 = 0 else n3 = nffeta(cc1,1/cc2,ier) endif if ( n3.ne.0 ) then print *,'ffdcrr: error: untested algorithm' ier = ier + 50 ipi12(1) = ipi12(1) + 4*12*n1*n3 endif endif else cc1 = cy*cfactz cc2 = cy*cfactw cs3(15) = (n1*zfflog(cc1,ieps1,czero,ier) + + n4*zfflog(cc2,ieps1,czero,ier))*c2ipi endif if ( n2.eq.n5 ) then if ( n2.eq.0 ) then * nothing to do else cc1 = cdwz/cdyz if ( absc(cc1) .lt. xloss ) then cs3(16) = n2*c2ipi*zfflo1(cc1,ier) else cc1 = -cdwy/cdyz cs3(16) = n2*c2ipi*zfflog(cc1,0,czero,ier) endif cc1 = -cy1*cfactz cc2 = -cy1*cfactw if ( Im(cc1).eq.0 .or. Im(cc2).eq.0 ) then n6 = 0 else n6 = nffeta(cc1,1/cc2,ier) endif if ( n6.ne.0 ) then print *,'ffdcrr: error: untested algorithm' ier = ier + 50 ipi12(2) = ipi12(2) + 4*12*n2*n6 endif endif else cc1 = -cy1*cfactz cc2 = -cy1*cfactw cs3(15) = (n2*zfflog(cc1,ieps2,czero,ier) + + n5*zfflog(cc2,ieps2,czero,ier))*c2ipi endif * #] only cancellations in cw, not in cy: * #[ Hill identity: elseif ( ( 1.gt.xloss*absc(cy) .or. absc(cc1).gt.xloss ) + .and. ( 1.gt.xloss*absc(cz) .or. absc(cz/cdyz).gt.xloss ) + .and. ( 1.gt.xloss*absc(cy) .or. absc(cdyz/cy).gt.xloss ) + ) then * do a Hill identity on the cy,cy-1 direction cyy = -cy*cw1/cdwy cyy1 = cw*cy1/cdwy czz = -cz*cw1/cdwz czz1 = cw*cz1/cdwz cdyyzz = -cw*cw1*(cdyz/(cdwy*cdwz)) ieps = -2*iepsz call ffcrr(cs3(1),ipi12(1),cyy,cyy1,czz,czz1,cdyyzz, + .FALSE.,czero,czero,czero,isoort,ieps,ier) cyy = cw1 cyy1 = cw czz = -cw1*cz/cdwz czz1 = cw*cz1/cdwz cdyyzz = cw*cw1/cdwz call ffcrr(cs3(8),ipi12(2),cyy,cyy1,czz,czz1,cdyyzz, + .FALSE.,czero,czero,czero,isoort,0,ier) do 30 i=8,14 - 30 cs3(i) = -cs3(i) + cs3(i) = -cs3(i) + 30 continue ipi12(2) = -ipi12(2) * the extra logarithms ... if ( 1 .lt. xloss*absc(cw) ) then chulp = zfflo1(1/cw,ier) else chulp = zfflog(-cw1/cw,0,czero,ier) endif cs3(15) = -zfflo1(cdwz/cdwy,ier)*chulp * #] Hill identity: * #[ Taylor expansion: else * Do a Taylor expansion if ( absc(cc1) .lt. xloss ) then cd3 = cdwz/cdwy * isign = 1 cc1n = cc1 cd2n = cd2 cd3n = cd3 cd21 = 1-cd2 cd21n1 = 1 cd31 = 1-cd3 cd31n1 = 1 ctot = cc1*cd2*cd3 do 50 i=2,20 cc1n = cc1n*cc1 cd21n1 = cd21n1*cd21 cd31n1 = cd31n1*cd31 cd2n = cd2n + cd2*cd21n1 cd3n = cd3n + cd3*cd31n1 cterm = cc1n*cd2n*cd3n*Re(xn2inv(i)) ctot = ctot + cterm if ( absc(cterm) .lt. precc*absc(ctot) ) goto 51 50 continue 51 continue cs3(1) = ctot elseif ( absc(cz/cdyz) .lt. xloss ) then call ffcrr(cs3(1),ipi12(1),cy,cy1,cz,cz1,cdyz,.TRUE., + cd2yzz,czp,czp1,isoort,iepsz,ier) call ffcrr(cs3(8),ipi12(2),cy,cy1,cw,cw1,-cdwy,.TRUE., + cd2yww,cwp,cwp1,isoort,iepsw,ier) do 110 i=8,14 - 110 cs3(i) = -cs3(i) + cs3(i) = -cs3(i) + 110 continue ipi12(2) = -ipi12(2) else call fferr(20,ier) return endif endif * #] Taylor expansion: *###] ffdcrr: end diff --git a/Shower/Dipole/Colorea/HelAmps_sm.cc b/Shower/Dipole/Colorea/HelAmps_sm.cc --- a/Shower/Dipole/Colorea/HelAmps_sm.cc +++ b/Shower/Dipole/Colorea/HelAmps_sm.cc @@ -1,1035 +1,1027 @@ //========================================================================== // This file has been automatically generated for C++ Standalone by // MadGraph5_aMC@NLO v. 2.5.4, 2017-03-28 // By the MadGraph5_aMC@NLO Development Team // Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch //========================================================================== #include "HelAmps_sm.h" #include #include #include #include using namespace std; namespace MG5_sm_COLOREA { void ixxxxx(double p[4], double fmass, int nhel, int nsf, complex fi[6]) { complex chi[2]; double sf[2], sfomega[2], omega[2], pp, pp3, sqp0p3, sqm[2]; int ip, im, nh; fi[0] = complex (-p[0] * nsf, -p[3] * nsf); fi[1] = complex (-p[1] * nsf, -p[2] * nsf); nh = nhel * nsf; if (fmass != 0.0) { pp = min(p[0], sqrt(p[1] * p[1] + p[2] * p[2] + p[3] * p[3])); if (pp == 0.0) { sqm[0] = sqrt(std::abs(fmass)); sqm[1] = Sgn(sqm[0], fmass); ip = (1 + nh)/2; im = (1 - nh)/2; fi[2] = ip * sqm[ip]; fi[3] = im * nsf * sqm[ip]; fi[4] = ip * nsf * sqm[im]; fi[5] = im * sqm[im]; } else { sf[0] = (1 + nsf + (1 - nsf) * nh) * 0.5; sf[1] = (1 + nsf - (1 - nsf) * nh) * 0.5; omega[0] = sqrt(p[0] + pp); omega[1] = fmass/omega[0]; ip = (1 + nh)/2; im = (1 - nh)/2; sfomega[0] = sf[0] * omega[ip]; sfomega[1] = sf[1] * omega[im]; pp3 = max(pp + p[3], 0.0); chi[0] = complex (sqrt(pp3 * 0.5/pp), 0); if (pp3 == 0.0) { chi[1] = complex (-nh, 0); } else { chi[1] = complex (nh * p[1], p[2])/sqrt(2.0 * pp * pp3); } fi[2] = sfomega[0] * chi[im]; fi[3] = sfomega[0] * chi[ip]; fi[4] = sfomega[1] * chi[im]; fi[5] = sfomega[1] * chi[ip]; } } else { if (p[1] == 0.0 and p[2] == 0.0 and p[3] < 0.0) { sqp0p3 = 0.0; } else { sqp0p3 = sqrt(max(p[0] + p[3], 0.0)) * nsf; } chi[0] = complex (sqp0p3, 0.0); if (sqp0p3 == 0.0) { chi[1] = complex (-nhel * sqrt(2.0 * p[0]), 0.0); } else { chi[1] = complex (nh * p[1], p[2])/sqp0p3; } if (nh == 1) { fi[2] = complex (0.0, 0.0); fi[3] = complex (0.0, 0.0); fi[4] = chi[0]; fi[5] = chi[1]; } else { fi[2] = chi[1]; fi[3] = chi[0]; fi[4] = complex (0.0, 0.0); fi[5] = complex (0.0, 0.0); } } return; } double Sgn(double a, double b) { return (b < 0)? - abs(a):abs(a); } void txxxxx(double p[4], double tmass, int nhel, int nst, complex tc[18]) { complex ft[6][4], ep[4], em[4], e0[4]; double pt, pt2, pp, pzpt, emp, sqh, sqs; int i, j; sqh = sqrt(0.5); sqs = sqrt(0.5/3); pt2 = p[1] * p[1] + p[2] * p[2]; pp = min(p[0], sqrt(pt2 + p[3] * p[3])); pt = min(pp, sqrt(pt2)); ft[4][0] = complex (p[0] * nst, p[3] * nst); ft[5][0] = complex (p[1] * nst, p[2] * nst); // construct eps+ if(nhel >= 0) { if(pp == 0) { ep[0] = complex (0, 0); ep[1] = complex (-sqh, 0); ep[2] = complex (0, nst * sqh); ep[3] = complex (0, 0); } else { ep[0] = complex (0, 0); ep[3] = complex (pt/pp * sqh, 0); if(pt != 0) { pzpt = p[3]/(pp * pt) * sqh; ep[1] = complex (-p[1] * pzpt, -nst * p[2]/pt * sqh); ep[2] = complex (-p[2] * pzpt, nst * p[1]/pt * sqh); } else { ep[1] = complex (-sqh, 0); ep[2] = complex (0, nst * Sgn(sqh, p[3])); } } } // construct eps- if(nhel <= 0) { if(pp == 0) { em[0] = complex (0, 0); em[1] = complex (sqh, 0); em[2] = complex (0, nst * sqh); em[3] = complex (0, 0); } else { em[0] = complex (0, 0); em[3] = complex (-pt/pp * sqh, 0); if(pt != 0) { pzpt = -p[3]/(pp * pt) * sqh; em[1] = complex (-p[1] * pzpt, -nst * p[2]/pt * sqh); em[2] = complex (-p[2] * pzpt, nst * p[1]/pt * sqh); } else { em[1] = complex (sqh, 0); em[2] = complex (0, nst * Sgn(sqh, p[3])); } } } // construct eps0 if(std::labs(nhel) <= 1) { if(pp == 0) { e0[0] = complex (0, 0); e0[1] = complex (0, 0); e0[2] = complex (0, 0); e0[3] = complex (1, 0); } else { emp = p[0]/(tmass * pp); e0[0] = complex (pp/tmass, 0); e0[3] = complex (p[3] * emp, 0); if(pt != 0) { e0[1] = complex (p[1] * emp, 0); e0[2] = complex (p[2] * emp, 0); } else { e0[1] = complex (0, 0); e0[2] = complex (0, 0); } } } if(nhel == 2) { for(j = 0; j < 4; j++ ) { for(i = 0; i < 4; i++ ) ft[i][j] = ep[i] * ep[j]; } } else if(nhel == -2) { for(j = 0; j < 4; j++ ) { for(i = 0; i < 4; i++ ) ft[i][j] = em[i] * em[j]; } } else if(tmass == 0) { for(j = 0; j < 4; j++ ) { for(i = 0; i < 4; i++ ) ft[i][j] = 0; } } else if(tmass != 0) { if(nhel == 1) { for(j = 0; j < 4; j++ ) { for(i = 0; i < 4; i++ ) ft[i][j] = sqh * (ep[i] * e0[j] + e0[i] * ep[j]); } } else if(nhel == 0) { for(j = 0; j < 4; j++ ) { for(i = 0; i < 4; i++ ) ft[i][j] = sqs * (ep[i] * em[j] + em[i] * ep[j] + 2.0 * e0[i] * e0[j]); } } else if(nhel == -1) { for(j = 0; j < 4; j++ ) { for(i = 0; i < 4; i++ ) ft[i][j] = sqh * (em[i] * e0[j] + e0[i] * em[j]); } } else { std::cerr << "Invalid helicity in txxxxx.\n"; std::exit(1); } } tc[0] = ft[4][0]; tc[1] = ft[5][0]; for(j = 0; j < 4; j++ ) { for(i = 0; i < 4; i++ ) tc[j * 4 + i + 2] = ft[j][i]; } } void vxxxxx(double p[4], double vmass, int nhel, int nsv, complex vc[6]) { double hel, hel0, pt, pt2, pp, pzpt, emp, sqh; int nsvahl; sqh = sqrt(0.5); hel = double(nhel); nsvahl = nsv * std::abs(hel); pt2 = (p[1] * p[1]) + (p[2] * p[2]); pp = min(p[0], sqrt(pt2 + (p[3] * p[3]))); pt = min(pp, sqrt(pt2)); vc[0] = complex (p[0] * nsv, p[3] * nsv); vc[1] = complex (p[1] * nsv, p[2] * nsv); if (vmass != 0.0) { hel0 = 1.0 - std::abs(hel); if(pp == 0.0) { vc[2] = complex (0.0, 0.0); vc[3] = complex (-hel * sqh, 0.0); vc[4] = complex (0.0, nsvahl * sqh); vc[5] = complex (hel0, 0.0); } else { emp = p[0]/(vmass * pp); vc[2] = complex (hel0 * pp/vmass, 0.0); vc[5] = complex (hel0 * p[3] * emp + hel * pt/pp * sqh, 0.0); if (pt != 0.0) { pzpt = p[3]/(pp * pt) * sqh * hel; vc[3] = complex (hel0 * p[1] * emp - p[1] * pzpt, -nsvahl * p[2]/pt * sqh); vc[4] = complex (hel0 * p[2] * emp - p[2] * pzpt, nsvahl * p[1]/pt * sqh); } else { vc[3] = complex (-hel * sqh, 0.0); vc[4] = complex (0.0, nsvahl * Sgn(sqh, p[3])); } } } else { pp = p[0]; pt = sqrt((p[1] * p[1]) + (p[2] * p[2])); vc[2] = complex (0.0, 0.0); vc[5] = complex (hel * pt/pp * sqh, 0.0); if (pt != 0.0) { pzpt = p[3]/(pp * pt) * sqh * hel; vc[3] = complex (-p[1] * pzpt, -nsv * p[2]/pt * sqh); vc[4] = complex (-p[2] * pzpt, nsv * p[1]/pt * sqh); } else { vc[3] = complex (-hel * sqh, 0.0); vc[4] = complex (0.0, nsv * Sgn(sqh, p[3])); } } return; } void sxxxxx(double p[4], int nss, complex sc[3]) { sc[2] = complex (1.00, 0.00); sc[0] = complex (p[0] * nss, p[3] * nss); sc[1] = complex (p[1] * nss, p[2] * nss); return; } void oxxxxx(double p[4], double fmass, int nhel, int nsf, complex fo[6]) { complex chi[2]; double sf[2], sfomeg[2], omega[2], pp, pp3, sqp0p3, sqm[2]; int nh, ip, im; fo[0] = complex (p[0] * nsf, p[3] * nsf); fo[1] = complex (p[1] * nsf, p[2] * nsf); nh = nhel * nsf; if (fmass != 0.000) { pp = min(p[0], sqrt((p[1] * p[1]) + (p[2] * p[2]) + (p[3] * p[3]))); if (pp == 0.000) { sqm[0] = sqrt(std::abs(fmass)); sqm[1] = Sgn(sqm[0], fmass); ip = -((1 - nh)/2) * nhel; im = (1 + nh)/2 * nhel; fo[2] = im * sqm[std::abs(ip)]; fo[3] = ip * nsf * sqm[std::abs(ip)]; fo[4] = im * nsf * sqm[std::abs(im)]; fo[5] = ip * sqm[std::abs(im)]; } else { pp = min(p[0], sqrt((p[1] * p[1]) + (p[2] * p[2]) + (p[3] * p[3]))); sf[0] = double(1 + nsf + (1 - nsf) * nh) * 0.5; sf[1] = double(1 + nsf - (1 - nsf) * nh) * 0.5; omega[0] = sqrt(p[0] + pp); omega[1] = fmass/omega[0]; ip = (1 + nh)/2; im = (1 - nh)/2; sfomeg[0] = sf[0] * omega[ip]; sfomeg[1] = sf[1] * omega[im]; pp3 = max(pp + p[3], 0.00); chi[0] = complex (sqrt(pp3 * 0.5/pp), 0.00); if (pp3 == 0.00) { chi[1] = complex (-nh, 0.00); } else { chi[1] = complex (nh * p[1], -p[2])/sqrt(2.0 * pp * pp3); } fo[2] = sfomeg[1] * chi[im]; fo[3] = sfomeg[1] * chi[ip]; fo[4] = sfomeg[0] * chi[im]; fo[5] = sfomeg[0] * chi[ip]; } } else { if((p[1] == 0.00) and (p[2] == 0.00) and (p[3] < 0.00)) { sqp0p3 = 0.00; } else { sqp0p3 = sqrt(max(p[0] + p[3], 0.00)) * nsf; } chi[0] = complex (sqp0p3, 0.00); if(sqp0p3 == 0.000) { chi[1] = complex (-nhel, 0.00) * sqrt(2.0 * p[0]); } else { chi[1] = complex (nh * p[1], -p[2])/sqp0p3; } if(nh == 1) { fo[2] = chi[0]; fo[3] = chi[1]; fo[4] = complex (0.00, 0.00); fo[5] = complex (0.00, 0.00); } else { fo[2] = complex (0.00, 0.00); fo[3] = complex (0.00, 0.00); fo[4] = chi[1]; fo[5] = chi[0]; } } return; } void FFV1P0_3(std::complex F1[], std::complex F2[], std::complex COUP, double M3, double W3, std::complex V3[]) { static std::complex cI = std::complex (0., 1.); double P3[4]; std::complex denom; V3[0] = +F1[0] + F2[0]; V3[1] = +F1[1] + F2[1]; P3[0] = -V3[0].real(); P3[1] = -V3[1].real(); P3[2] = -V3[1].imag(); P3[3] = -V3[0].imag(); denom = COUP/((P3[0] * P3[0]) - (P3[1] * P3[1]) - (P3[2] * P3[2]) - (P3[3] * P3[3]) - M3 * (M3 - cI * W3)); V3[2] = denom * (-cI) * (F1[2] * F2[4] + F1[3] * F2[5] + F1[4] * F2[2] + F1[5] * F2[3]); V3[3] = denom * (-cI) * (F1[4] * F2[3] + F1[5] * F2[2] - F1[2] * F2[5] - F1[3] * F2[4]); V3[4] = denom * (-cI) * (-cI * (F1[2] * F2[5] + F1[5] * F2[2]) + cI * (F1[3] * F2[4] + F1[4] * F2[3])); V3[5] = denom * (-cI) * (F1[3] * F2[5] + F1[4] * F2[2] - F1[2] * F2[4] - F1[5] * F2[3]); } void FFV2_2(std::complex F1[], std::complex V3[], std::complex COUP, double M2, double W2, std::complex F2[]) { static std::complex cI = std::complex (0., 1.); double P2[4]; std::complex denom; F2[0] = +F1[0] + V3[0]; F2[1] = +F1[1] + V3[1]; P2[0] = -F2[0].real(); P2[1] = -F2[1].real(); P2[2] = -F2[1].imag(); P2[3] = -F2[0].imag(); denom = COUP/((P2[0] * P2[0]) - (P2[1] * P2[1]) - (P2[2] * P2[2]) - (P2[3] * P2[3]) - M2 * (M2 - cI * W2)); F2[2] = denom * cI * (F1[2] * (P2[0] * (V3[2] + V3[5]) + (P2[1] * (-1.) * (V3[3] + cI * (V3[4])) + (P2[2] * (+cI * (V3[3]) - V3[4]) - P2[3] * (V3[2] + V3[5])))) + F1[3] * (P2[0] * (V3[3] - cI * (V3[4])) + (P2[1] * (V3[5] - V3[2]) + (P2[2] * (-cI * (V3[5]) + cI * (V3[2])) + P2[3] * (+cI * (V3[4]) - V3[3]))))); F2[3] = denom * cI * (F1[2] * (P2[0] * (V3[3] + cI * (V3[4])) + (P2[1] * (-1.) * (V3[2] + V3[5]) + (P2[2] * (-1.) * (+cI * (V3[2] + V3[5])) + P2[3] * (V3[3] + cI * (V3[4]))))) + F1[3] * (P2[0] * (V3[2] - V3[5]) + (P2[1] * (+cI * (V3[4]) - V3[3]) + (P2[2] * (-1.) * (V3[4] + cI * (V3[3])) + P2[3] * (V3[2] - V3[5]))))); F2[4] = denom * - cI * M2 * (F1[2] * (-1.) * (V3[2] + V3[5]) + F1[3] * (+cI * (V3[4]) - V3[3])); F2[5] = denom * cI * M2 * (F1[2] * (V3[3] + cI * (V3[4])) + F1[3] * (V3[2] - V3[5])); } void FFV2_5_2(std::complex F1[], std::complex V3[], std::complex COUP1, std::complex COUP2, double M2, double W2, std::complex F2[]) { - static std::complex cI = std::complex (0., 1.); - std::complex Ftmp[6]; - double P2[4]; + std::complex Ftmp[6]; std::complex denom; int i; FFV2_2(F1, V3, COUP1, M2, W2, F2); FFV5_2(F1, V3, COUP2, M2, W2, Ftmp); i = 2; while (i < 6) { F2[i] = F2[i] + Ftmp[i]; i++; } } void FFV1_2(std::complex F1[], std::complex V3[], std::complex COUP, double M2, double W2, std::complex F2[]) { static std::complex cI = std::complex (0., 1.); double P2[4]; std::complex denom; F2[0] = +F1[0] + V3[0]; F2[1] = +F1[1] + V3[1]; P2[0] = -F2[0].real(); P2[1] = -F2[1].real(); P2[2] = -F2[1].imag(); P2[3] = -F2[0].imag(); denom = COUP/((P2[0] * P2[0]) - (P2[1] * P2[1]) - (P2[2] * P2[2]) - (P2[3] * P2[3]) - M2 * (M2 - cI * W2)); F2[2] = denom * cI * (F1[2] * (P2[0] * (V3[2] + V3[5]) + (P2[1] * (-1.) * (V3[3] + cI * (V3[4])) + (P2[2] * (+cI * (V3[3]) - V3[4]) - P2[3] * (V3[2] + V3[5])))) + (F1[3] * (P2[0] * (V3[3] - cI * (V3[4])) + (P2[1] * (V3[5] - V3[2]) + (P2[2] * (-cI * (V3[5]) + cI * (V3[2])) + P2[3] * (+cI * (V3[4]) - V3[3])))) + M2 * (F1[4] * (V3[2] - V3[5]) + F1[5] * (+cI * (V3[4]) - V3[3])))); F2[3] = denom * (-cI) * (F1[2] * (P2[0] * (-1.) * (V3[3] + cI * (V3[4])) + (P2[1] * (V3[2] + V3[5]) + (P2[2] * (+cI * (V3[2] + V3[5])) - P2[3] * (V3[3] + cI * (V3[4]))))) + (F1[3] * (P2[0] * (V3[5] - V3[2]) + (P2[1] * (V3[3] - cI * (V3[4])) + (P2[2] * (V3[4] + cI * (V3[3])) + P2[3] * (V3[5] - V3[2])))) + M2 * (F1[4] * (V3[3] + cI * (V3[4])) - F1[5] * (V3[2] + V3[5])))); F2[4] = denom * (-cI) * (F1[4] * (P2[0] * (V3[5] - V3[2]) + (P2[1] * (V3[3] + cI * (V3[4])) + (P2[2] * (V3[4] - cI * (V3[3])) + P2[3] * (V3[5] - V3[2])))) + (F1[5] * (P2[0] * (V3[3] - cI * (V3[4])) + (P2[1] * (-1.) * (V3[2] + V3[5]) + (P2[2] * (+cI * (V3[2] + V3[5])) + P2[3] * (V3[3] - cI * (V3[4]))))) + M2 * (F1[2] * (-1.) * (V3[2] + V3[5]) + F1[3] * (+cI * (V3[4]) - V3[3])))); F2[5] = denom * cI * (F1[4] * (P2[0] * (-1.) * (V3[3] + cI * (V3[4])) + (P2[1] * (V3[2] - V3[5]) + (P2[2] * (-cI * (V3[5]) + cI * (V3[2])) + P2[3] * (V3[3] + cI * (V3[4]))))) + (F1[5] * (P2[0] * (V3[2] + V3[5]) + (P2[1] * (+cI * (V3[4]) - V3[3]) + (P2[2] * (-1.) * (V3[4] + cI * (V3[3])) - P2[3] * (V3[2] + V3[5])))) + M2 * (F1[2] * (V3[3] + cI * (V3[4])) + F1[3] * (V3[2] - V3[5])))); } void FFV2_0(std::complex F1[], std::complex F2[], std::complex V3[], std::complex COUP, std::complex & vertex) { static std::complex cI = std::complex (0., 1.); std::complex TMP9; TMP9 = (F1[2] * (F2[4] * (V3[2] + V3[5]) + F2[5] * (V3[3] + cI * (V3[4]))) + F1[3] * (F2[4] * (V3[3] - cI * (V3[4])) + F2[5] * (V3[2] - V3[5]))); vertex = COUP * - cI * TMP9; } void FFV2_5_0(std::complex F1[], std::complex F2[], std::complex V3[], std::complex COUP1, std::complex COUP2, std::complex & vertex) { - static std::complex cI = std::complex (0., 1.); std::complex tmp; FFV2_0(F1, F2, V3, COUP1, vertex); FFV5_0(F1, F2, V3, COUP2, tmp); vertex = vertex + tmp; } void FFV5_1(std::complex F2[], std::complex V3[], std::complex COUP, double M1, double W1, std::complex F1[]) { static std::complex cI = std::complex (0., 1.); double P1[4]; std::complex denom; F1[0] = +F2[0] + V3[0]; F1[1] = +F2[1] + V3[1]; P1[0] = -F1[0].real(); P1[1] = -F1[1].real(); P1[2] = -F1[1].imag(); P1[3] = -F1[0].imag(); denom = COUP/((P1[0] * P1[0]) - (P1[1] * P1[1]) - (P1[2] * P1[2]) - (P1[3] * P1[3]) - M1 * (M1 - cI * W1)); F1[2] = denom * 4. * cI * (F2[2] * (P1[0] * (V3[5] - V3[2]) + (P1[1] * (V3[3] - cI * (V3[4])) + (P1[2] * (V3[4] + cI * (V3[3])) + P1[3] * (V3[5] - V3[2])))) + (+1./4. * (M1 * (F2[5] * (V3[3] + cI * (V3[4])) + 4. * (F2[4] * 1./4. * (V3[2] + V3[5])))) + F2[3] * (P1[0] * (V3[3] + cI * (V3[4])) + (P1[1] * (-1.) * (V3[2] + V3[5]) + (P1[2] * (-1.) * (+cI * (V3[2] + V3[5])) + P1[3] * (V3[3] + cI * (V3[4]))))))); F1[3] = denom * 4. * cI * (F2[2] * (P1[0] * (V3[3] - cI * (V3[4])) + (P1[1] * (V3[5] - V3[2]) + (P1[2] * (-cI * (V3[5]) + cI * (V3[2])) + P1[3] * (+cI * (V3[4]) - V3[3])))) + (+1./4. * (M1 * (F2[5] * (V3[2] - V3[5]) + 4. * (F2[4] * 1./4. * (V3[3] - cI * (V3[4]))))) + F2[3] * (P1[0] * (-1.) * (V3[2] + V3[5]) + (P1[1] * (V3[3] + cI * (V3[4])) + (P1[2] * (V3[4] - cI * (V3[3])) + P1[3] * (V3[2] + V3[5])))))); F1[4] = denom * (-cI) * (F2[4] * (P1[0] * (V3[2] + V3[5]) + (P1[1] * (+cI * (V3[4]) - V3[3]) + (P1[2] * (-1.) * (V3[4] + cI * (V3[3])) - P1[3] * (V3[2] + V3[5])))) + (F2[5] * (P1[0] * (V3[3] + cI * (V3[4])) + (P1[1] * (V3[5] - V3[2]) + (P1[2] * (-cI * (V3[2]) + cI * (V3[5])) - P1[3] * (V3[3] + cI * (V3[4]))))) + M1 * (F2[2] * 4. * (V3[5] - V3[2]) + 4. * (F2[3] * (V3[3] + cI * (V3[4])))))); F1[5] = denom * cI * (F2[4] * (P1[0] * (+cI * (V3[4]) - V3[3]) + (P1[1] * (V3[2] + V3[5]) + (P1[2] * (-1.) * (+cI * (V3[2] + V3[5])) + P1[3] * (+cI * (V3[4]) - V3[3])))) + (F2[5] * (P1[0] * (V3[5] - V3[2]) + (P1[1] * (V3[3] + cI * (V3[4])) + (P1[2] * (V3[4] - cI * (V3[3])) + P1[3] * (V3[5] - V3[2])))) + M1 * (F2[2] * 4. * (+cI * (V3[4]) - V3[3]) + 4. * (F2[3] * (V3[2] + V3[5]))))); } void FFV1_0(std::complex F1[], std::complex F2[], std::complex V3[], std::complex COUP, std::complex & vertex) { static std::complex cI = std::complex (0., 1.); std::complex TMP13; TMP13 = (F1[2] * (F2[4] * (V3[2] + V3[5]) + F2[5] * (V3[3] + cI * (V3[4]))) + (F1[3] * (F2[4] * (V3[3] - cI * (V3[4])) + F2[5] * (V3[2] - V3[5])) + (F1[4] * (F2[2] * (V3[2] - V3[5]) - F2[3] * (V3[3] + cI * (V3[4]))) + F1[5] * (F2[2] * (+cI * (V3[4]) - V3[3]) + F2[3] * (V3[2] + V3[5]))))); vertex = COUP * - cI * TMP13; } void VVVV4P0_1(std::complex V2[], std::complex V3[], std::complex V4[], std::complex COUP, double M1, double W1, std::complex V1[]) { static std::complex cI = std::complex (0., 1.); std::complex TMP12; std::complex TMP11; double P1[4]; std::complex denom; V1[0] = +V2[0] + V3[0] + V4[0]; V1[1] = +V2[1] + V3[1] + V4[1]; P1[0] = -V1[0].real(); P1[1] = -V1[1].real(); P1[2] = -V1[1].imag(); P1[3] = -V1[0].imag(); TMP11 = (V2[2] * V4[2] - V2[3] * V4[3] - V2[4] * V4[4] - V2[5] * V4[5]); TMP12 = (V3[2] * V4[2] - V3[3] * V4[3] - V3[4] * V4[4] - V3[5] * V4[5]); denom = COUP/((P1[0] * P1[0]) - (P1[1] * P1[1]) - (P1[2] * P1[2]) - (P1[3] * P1[3]) - M1 * (M1 - cI * W1)); V1[2] = denom * (-cI * (V3[2] * TMP11) + cI * (V2[2] * TMP12)); V1[3] = denom * (-cI * (V3[3] * TMP11) + cI * (V2[3] * TMP12)); V1[4] = denom * (-cI * (V3[4] * TMP11) + cI * (V2[4] * TMP12)); V1[5] = denom * (-cI * (V3[5] * TMP11) + cI * (V2[5] * TMP12)); } void VVVV3P0_1(std::complex V2[], std::complex V3[], std::complex V4[], std::complex COUP, double M1, double W1, std::complex V1[]) { static std::complex cI = std::complex (0., 1.); std::complex TMP12; double P1[4]; std::complex TMP6; std::complex denom; V1[0] = +V2[0] + V3[0] + V4[0]; V1[1] = +V2[1] + V3[1] + V4[1]; P1[0] = -V1[0].real(); P1[1] = -V1[1].real(); P1[2] = -V1[1].imag(); P1[3] = -V1[0].imag(); TMP6 = (V3[2] * V2[2] - V3[3] * V2[3] - V3[4] * V2[4] - V3[5] * V2[5]); TMP12 = (V3[2] * V4[2] - V3[3] * V4[3] - V3[4] * V4[4] - V3[5] * V4[5]); denom = COUP/((P1[0] * P1[0]) - (P1[1] * P1[1]) - (P1[2] * P1[2]) - (P1[3] * P1[3]) - M1 * (M1 - cI * W1)); V1[2] = denom * (-cI * (TMP6 * V4[2]) + cI * (V2[2] * TMP12)); V1[3] = denom * (-cI * (TMP6 * V4[3]) + cI * (V2[3] * TMP12)); V1[4] = denom * (-cI * (TMP6 * V4[4]) + cI * (V2[4] * TMP12)); V1[5] = denom * (-cI * (TMP6 * V4[5]) + cI * (V2[5] * TMP12)); } void VVV1_0(std::complex V1[], std::complex V2[], std::complex V3[], std::complex COUP, std::complex & vertex) { static std::complex cI = std::complex (0., 1.); std::complex TMP2; std::complex TMP1; double P1[4]; std::complex TMP0; double P2[4]; std::complex TMP7; double P3[4]; std::complex TMP6; std::complex TMP5; std::complex TMP4; std::complex TMP3; std::complex TMP8; P1[0] = V1[0].real(); P1[1] = V1[1].real(); P1[2] = V1[1].imag(); P1[3] = V1[0].imag(); P2[0] = V2[0].real(); P2[1] = V2[1].real(); P2[2] = V2[1].imag(); P2[3] = V2[0].imag(); P3[0] = V3[0].real(); P3[1] = V3[1].real(); P3[2] = V3[1].imag(); P3[3] = V3[0].imag(); TMP8 = (V1[2] * P3[0] - V1[3] * P3[1] - V1[4] * P3[2] - V1[5] * P3[3]); TMP5 = (V2[2] * P3[0] - V2[3] * P3[1] - V2[4] * P3[2] - V2[5] * P3[3]); TMP4 = (P1[0] * V2[2] - P1[1] * V2[3] - P1[2] * V2[4] - P1[3] * V2[5]); TMP7 = (V1[2] * P2[0] - V1[3] * P2[1] - V1[4] * P2[2] - V1[5] * P2[3]); TMP6 = (V3[2] * V2[2] - V3[3] * V2[3] - V3[4] * V2[4] - V3[5] * V2[5]); TMP1 = (V2[2] * V1[2] - V2[3] * V1[3] - V2[4] * V1[4] - V2[5] * V1[5]); TMP0 = (V3[2] * P1[0] - V3[3] * P1[1] - V3[4] * P1[2] - V3[5] * P1[3]); TMP3 = (V3[2] * V1[2] - V3[3] * V1[3] - V3[4] * V1[4] - V3[5] * V1[5]); TMP2 = (V3[2] * P2[0] - V3[3] * P2[1] - V3[4] * P2[2] - V3[5] * P2[3]); vertex = COUP * (TMP1 * (-cI * (TMP0) + cI * (TMP2)) + (TMP3 * (-cI * (TMP5) + cI * (TMP4)) + TMP6 * (-cI * (TMP7) + cI * (TMP8)))); } void FFV2_3(std::complex F1[], std::complex F2[], std::complex COUP, double M3, double W3, std::complex V3[]) { static std::complex cI = std::complex (0., 1.); std::complex denom; std::complex TMP10; double P3[4]; double OM3; OM3 = 0.; if (M3 != 0.) OM3 = 1./(M3 * M3); V3[0] = +F1[0] + F2[0]; V3[1] = +F1[1] + F2[1]; P3[0] = -V3[0].real(); P3[1] = -V3[1].real(); P3[2] = -V3[1].imag(); P3[3] = -V3[0].imag(); TMP10 = (F1[2] * (F2[4] * (P3[0] + P3[3]) + F2[5] * (P3[1] + cI * (P3[2]))) + F1[3] * (F2[4] * (P3[1] - cI * (P3[2])) + F2[5] * (P3[0] - P3[3]))); denom = COUP/((P3[0] * P3[0]) - (P3[1] * P3[1]) - (P3[2] * P3[2]) - (P3[3] * P3[3]) - M3 * (M3 - cI * W3)); V3[2] = denom * (-cI) * (F1[2] * F2[4] + F1[3] * F2[5] - P3[0] * OM3 * TMP10); V3[3] = denom * (-cI) * (-F1[2] * F2[5] - F1[3] * F2[4] - P3[1] * OM3 * TMP10); V3[4] = denom * (-cI) * (-cI * (F1[2] * F2[5]) + cI * (F1[3] * F2[4]) - P3[2] * OM3 * TMP10); V3[5] = denom * (-cI) * (F1[3] * F2[5] - F1[2] * F2[4] - P3[3] * OM3 * TMP10); } void FFV2_4_3(std::complex F1[], std::complex F2[], std::complex COUP1, std::complex COUP2, double M3, double W3, std::complex V3[]) { - static std::complex cI = std::complex (0., 1.); - std::complex denom; - double P3[4]; - double OM3; + std::complex denom; int i; std::complex Vtmp[6]; FFV2_3(F1, F2, COUP1, M3, W3, V3); FFV4_3(F1, F2, COUP2, M3, W3, Vtmp); i = 2; while (i < 6) { V3[i] = V3[i] + Vtmp[i]; i++; } } void FFV5_2(std::complex F1[], std::complex V3[], std::complex COUP, double M2, double W2, std::complex F2[]) { static std::complex cI = std::complex (0., 1.); double P2[4]; std::complex denom; F2[0] = +F1[0] + V3[0]; F2[1] = +F1[1] + V3[1]; P2[0] = -F2[0].real(); P2[1] = -F2[1].real(); P2[2] = -F2[1].imag(); P2[3] = -F2[0].imag(); denom = COUP/((P2[0] * P2[0]) - (P2[1] * P2[1]) - (P2[2] * P2[2]) - (P2[3] * P2[3]) - M2 * (M2 - cI * W2)); F2[2] = denom * cI * (F1[2] * (P2[0] * (V3[2] + V3[5]) + (P2[1] * (-1.) * (V3[3] + cI * (V3[4])) + (P2[2] * (+cI * (V3[3]) - V3[4]) - P2[3] * (V3[2] + V3[5])))) + (F1[3] * (P2[0] * (V3[3] - cI * (V3[4])) + (P2[1] * (V3[5] - V3[2]) + (P2[2] * (-cI * (V3[5]) + cI * (V3[2])) + P2[3] * (+cI * (V3[4]) - V3[3])))) + M2 * (F1[4] * 4. * (V3[2] - V3[5]) + 4. * (F1[5] * (+cI * (V3[4]) - V3[3]))))); F2[3] = denom * cI * (F1[2] * (P2[0] * (V3[3] + cI * (V3[4])) + (P2[1] * (-1.) * (V3[2] + V3[5]) + (P2[2] * (-1.) * (+cI * (V3[2] + V3[5])) + P2[3] * (V3[3] + cI * (V3[4]))))) + (F1[3] * (P2[0] * (V3[2] - V3[5]) + (P2[1] * (+cI * (V3[4]) - V3[3]) + (P2[2] * (-1.) * (V3[4] + cI * (V3[3])) + P2[3] * (V3[2] - V3[5])))) + M2 * (F1[4] * (-4.) * (V3[3] + cI * (V3[4])) + 4. * (F1[5] * (V3[2] + V3[5]))))); F2[4] = denom * (-4. * cI) * (F1[4] * (P2[0] * (V3[5] - V3[2]) + (P2[1] * (V3[3] + cI * (V3[4])) + (P2[2] * (V3[4] - cI * (V3[3])) + P2[3] * (V3[5] - V3[2])))) + (+1./4. * (M2 * (F1[3] * (+cI * (V3[4]) - V3[3]) + 4. * (F1[2] * (-1./4.) * (V3[2] + V3[5])))) + F1[5] * (P2[0] * (V3[3] - cI * (V3[4])) + (P2[1] * (-1.) * (V3[2] + V3[5]) + (P2[2] * (+cI * (V3[2] + V3[5])) + P2[3] * (V3[3] - cI * (V3[4]))))))); F2[5] = denom * (-4. * cI) * (F1[4] * (P2[0] * (V3[3] + cI * (V3[4])) + (P2[1] * (V3[5] - V3[2]) + (P2[2] * (-cI * (V3[2]) + cI * (V3[5])) - P2[3] * (V3[3] + cI * (V3[4]))))) + (+1./4. * (M2 * (F1[3] * (V3[5] - V3[2]) + 4. * (F1[2] * (-1./4.) * (V3[3] + cI * (V3[4]))))) + F1[5] * (P2[0] * (-1.) * (V3[2] + V3[5]) + (P2[1] * (V3[3] - cI * (V3[4])) + (P2[2] * (V3[4] + cI * (V3[3])) + P2[3] * (V3[2] + V3[5])))))); } void FFV2_1(std::complex F2[], std::complex V3[], std::complex COUP, double M1, double W1, std::complex F1[]) { static std::complex cI = std::complex (0., 1.); double P1[4]; std::complex denom; F1[0] = +F2[0] + V3[0]; F1[1] = +F2[1] + V3[1]; P1[0] = -F1[0].real(); P1[1] = -F1[1].real(); P1[2] = -F1[1].imag(); P1[3] = -F1[0].imag(); denom = COUP/((P1[0] * P1[0]) - (P1[1] * P1[1]) - (P1[2] * P1[2]) - (P1[3] * P1[3]) - M1 * (M1 - cI * W1)); F1[2] = denom * cI * M1 * (F2[4] * (V3[2] + V3[5]) + F2[5] * (V3[3] + cI * (V3[4]))); F1[3] = denom * - cI * M1 * (F2[4] * (+cI * (V3[4]) - V3[3]) + F2[5] * (V3[5] - V3[2])); F1[4] = denom * (-cI) * (F2[4] * (P1[0] * (V3[2] + V3[5]) + (P1[1] * (+cI * (V3[4]) - V3[3]) + (P1[2] * (-1.) * (V3[4] + cI * (V3[3])) - P1[3] * (V3[2] + V3[5])))) + F2[5] * (P1[0] * (V3[3] + cI * (V3[4])) + (P1[1] * (V3[5] - V3[2]) + (P1[2] * (-cI * (V3[2]) + cI * (V3[5])) - P1[3] * (V3[3] + cI * (V3[4])))))); F1[5] = denom * (-cI) * (F2[4] * (P1[0] * (V3[3] - cI * (V3[4])) + (P1[1] * (-1.) * (V3[2] + V3[5]) + (P1[2] * (+cI * (V3[2] + V3[5])) + P1[3] * (V3[3] - cI * (V3[4]))))) + F2[5] * (P1[0] * (V3[2] - V3[5]) + (P1[1] * (-1.) * (V3[3] + cI * (V3[4])) + (P1[2] * (+cI * (V3[3]) - V3[4]) + P1[3] * (V3[2] - V3[5]))))); } void FFV2_5_1(std::complex F2[], std::complex V3[], std::complex COUP1, std::complex COUP2, double M1, double W1, std::complex F1[]) { - static std::complex cI = std::complex (0., 1.); - double P1[4]; std::complex denom; int i; std::complex Ftmp[6]; FFV2_1(F2, V3, COUP1, M1, W1, F1); FFV5_1(F2, V3, COUP2, M1, W1, Ftmp); i = 2; while (i < 6) { F1[i] = F1[i] + Ftmp[i]; i++; } } void FFV5_0(std::complex F1[], std::complex F2[], std::complex V3[], std::complex COUP, std::complex & vertex) { static std::complex cI = std::complex (0., 1.); std::complex TMP15; std::complex TMP16; TMP15 = (F1[2] * (F2[4] * (V3[2] + V3[5]) + F2[5] * (V3[3] + cI * (V3[4]))) + F1[3] * (F2[4] * (V3[3] - cI * (V3[4])) + F2[5] * (V3[2] - V3[5]))); TMP16 = (F1[4] * (F2[2] * (V3[2] - V3[5]) - F2[3] * (V3[3] + cI * (V3[4]))) + F1[5] * (F2[2] * (+cI * (V3[4]) - V3[3]) + F2[3] * (V3[2] + V3[5]))); vertex = COUP * (-1.) * (+cI * (TMP15) + 4. * cI * (TMP16)); } void FFV1_1(std::complex F2[], std::complex V3[], std::complex COUP, double M1, double W1, std::complex F1[]) { static std::complex cI = std::complex (0., 1.); double P1[4]; std::complex denom; F1[0] = +F2[0] + V3[0]; F1[1] = +F2[1] + V3[1]; P1[0] = -F1[0].real(); P1[1] = -F1[1].real(); P1[2] = -F1[1].imag(); P1[3] = -F1[0].imag(); denom = COUP/((P1[0] * P1[0]) - (P1[1] * P1[1]) - (P1[2] * P1[2]) - (P1[3] * P1[3]) - M1 * (M1 - cI * W1)); F1[2] = denom * cI * (F2[2] * (P1[0] * (V3[5] - V3[2]) + (P1[1] * (V3[3] - cI * (V3[4])) + (P1[2] * (V3[4] + cI * (V3[3])) + P1[3] * (V3[5] - V3[2])))) + (F2[3] * (P1[0] * (V3[3] + cI * (V3[4])) + (P1[1] * (-1.) * (V3[2] + V3[5]) + (P1[2] * (-1.) * (+cI * (V3[2] + V3[5])) + P1[3] * (V3[3] + cI * (V3[4]))))) + M1 * (F2[4] * (V3[2] + V3[5]) + F2[5] * (V3[3] + cI * (V3[4]))))); F1[3] = denom * (-cI) * (F2[2] * (P1[0] * (+cI * (V3[4]) - V3[3]) + (P1[1] * (V3[2] - V3[5]) + (P1[2] * (-cI * (V3[2]) + cI * (V3[5])) + P1[3] * (V3[3] - cI * (V3[4]))))) + (F2[3] * (P1[0] * (V3[2] + V3[5]) + (P1[1] * (-1.) * (V3[3] + cI * (V3[4])) + (P1[2] * (+cI * (V3[3]) - V3[4]) - P1[3] * (V3[2] + V3[5])))) + M1 * (F2[4] * (+cI * (V3[4]) - V3[3]) + F2[5] * (V3[5] - V3[2])))); F1[4] = denom * (-cI) * (F2[4] * (P1[0] * (V3[2] + V3[5]) + (P1[1] * (+cI * (V3[4]) - V3[3]) + (P1[2] * (-1.) * (V3[4] + cI * (V3[3])) - P1[3] * (V3[2] + V3[5])))) + (F2[5] * (P1[0] * (V3[3] + cI * (V3[4])) + (P1[1] * (V3[5] - V3[2]) + (P1[2] * (-cI * (V3[2]) + cI * (V3[5])) - P1[3] * (V3[3] + cI * (V3[4]))))) + M1 * (F2[2] * (V3[5] - V3[2]) + F2[3] * (V3[3] + cI * (V3[4]))))); F1[5] = denom * cI * (F2[4] * (P1[0] * (+cI * (V3[4]) - V3[3]) + (P1[1] * (V3[2] + V3[5]) + (P1[2] * (-1.) * (+cI * (V3[2] + V3[5])) + P1[3] * (+cI * (V3[4]) - V3[3])))) + (F2[5] * (P1[0] * (V3[5] - V3[2]) + (P1[1] * (V3[3] + cI * (V3[4])) + (P1[2] * (V3[4] - cI * (V3[3])) + P1[3] * (V3[5] - V3[2])))) + M1 * (F2[2] * (+cI * (V3[4]) - V3[3]) + F2[3] * (V3[2] + V3[5])))); } void FFV4_3(std::complex F1[], std::complex F2[], std::complex COUP, double M3, double W3, std::complex V3[]) { static std::complex cI = std::complex (0., 1.); std::complex denom; std::complex TMP10; double P3[4]; double OM3; std::complex TMP14; OM3 = 0.; if (M3 != 0.) OM3 = 1./(M3 * M3); V3[0] = +F1[0] + F2[0]; V3[1] = +F1[1] + F2[1]; P3[0] = -V3[0].real(); P3[1] = -V3[1].real(); P3[2] = -V3[1].imag(); P3[3] = -V3[0].imag(); TMP14 = (F1[4] * (F2[2] * (P3[0] - P3[3]) - F2[3] * (P3[1] + cI * (P3[2]))) + F1[5] * (F2[2] * (+cI * (P3[2]) - P3[1]) + F2[3] * (P3[0] + P3[3]))); TMP10 = (F1[2] * (F2[4] * (P3[0] + P3[3]) + F2[5] * (P3[1] + cI * (P3[2]))) + F1[3] * (F2[4] * (P3[1] - cI * (P3[2])) + F2[5] * (P3[0] - P3[3]))); denom = COUP/((P3[0] * P3[0]) - (P3[1] * P3[1]) - (P3[2] * P3[2]) - (P3[3] * P3[3]) - M3 * (M3 - cI * W3)); V3[2] = denom * (-2. * cI) * (OM3 * - 1./2. * P3[0] * (TMP10 + 2. * (TMP14)) + (+1./2. * (F1[2] * F2[4] + F1[3] * F2[5]) + F1[4] * F2[2] + F1[5] * F2[3])); V3[3] = denom * (-2. * cI) * (OM3 * - 1./2. * P3[1] * (TMP10 + 2. * (TMP14)) + (-1./2. * (F1[2] * F2[5] + F1[3] * F2[4]) + F1[4] * F2[3] + F1[5] * F2[2])); V3[4] = denom * 2. * cI * (OM3 * 1./2. * P3[2] * (TMP10 + 2. * (TMP14)) + (+1./2. * cI * (F1[2] * F2[5]) - 1./2. * cI * (F1[3] * F2[4]) - cI * (F1[4] * F2[3]) + cI * (F1[5] * F2[2]))); V3[5] = denom * 2. * cI * (OM3 * 1./2. * P3[3] * (TMP10 + 2. * (TMP14)) + (+1./2. * (F1[2] * F2[4]) - 1./2. * (F1[3] * F2[5]) - F1[4] * F2[2] + F1[5] * F2[3])); } void VVVV1P0_1(std::complex V2[], std::complex V3[], std::complex V4[], std::complex COUP, double M1, double W1, std::complex V1[]) { static std::complex cI = std::complex (0., 1.); std::complex TMP11; double P1[4]; std::complex TMP6; std::complex denom; V1[0] = +V2[0] + V3[0] + V4[0]; V1[1] = +V2[1] + V3[1] + V4[1]; P1[0] = -V1[0].real(); P1[1] = -V1[1].real(); P1[2] = -V1[1].imag(); P1[3] = -V1[0].imag(); TMP6 = (V3[2] * V2[2] - V3[3] * V2[3] - V3[4] * V2[4] - V3[5] * V2[5]); TMP11 = (V2[2] * V4[2] - V2[3] * V4[3] - V2[4] * V4[4] - V2[5] * V4[5]); denom = COUP/((P1[0] * P1[0]) - (P1[1] * P1[1]) - (P1[2] * P1[2]) - (P1[3] * P1[3]) - M1 * (M1 - cI * W1)); V1[2] = denom * (-cI * (TMP6 * V4[2]) + cI * (V3[2] * TMP11)); V1[3] = denom * (-cI * (TMP6 * V4[3]) + cI * (V3[3] * TMP11)); V1[4] = denom * (-cI * (TMP6 * V4[4]) + cI * (V3[4] * TMP11)); V1[5] = denom * (-cI * (TMP6 * V4[5]) + cI * (V3[5] * TMP11)); } void VVV1P0_1(std::complex V2[], std::complex V3[], std::complex COUP, double M1, double W1, std::complex V1[]) { static std::complex cI = std::complex (0., 1.); std::complex TMP2; double P1[4]; std::complex TMP0; double P2[4]; double P3[4]; std::complex TMP6; std::complex TMP5; std::complex TMP4; std::complex denom; P2[0] = V2[0].real(); P2[1] = V2[1].real(); P2[2] = V2[1].imag(); P2[3] = V2[0].imag(); P3[0] = V3[0].real(); P3[1] = V3[1].real(); P3[2] = V3[1].imag(); P3[3] = V3[0].imag(); V1[0] = +V2[0] + V3[0]; V1[1] = +V2[1] + V3[1]; P1[0] = -V1[0].real(); P1[1] = -V1[1].real(); P1[2] = -V1[1].imag(); P1[3] = -V1[0].imag(); TMP5 = (V2[2] * P3[0] - V2[3] * P3[1] - V2[4] * P3[2] - V2[5] * P3[3]); TMP4 = (P1[0] * V2[2] - P1[1] * V2[3] - P1[2] * V2[4] - P1[3] * V2[5]); TMP6 = (V3[2] * V2[2] - V3[3] * V2[3] - V3[4] * V2[4] - V3[5] * V2[5]); TMP0 = (V3[2] * P1[0] - V3[3] * P1[1] - V3[4] * P1[2] - V3[5] * P1[3]); TMP2 = (V3[2] * P2[0] - V3[3] * P2[1] - V3[4] * P2[2] - V3[5] * P2[3]); denom = COUP/((P1[0] * P1[0]) - (P1[1] * P1[1]) - (P1[2] * P1[2]) - (P1[3] * P1[3]) - M1 * (M1 - cI * W1)); V1[2] = denom * (TMP6 * (-cI * (P2[0]) + cI * (P3[0])) + (V2[2] * (-cI * (TMP0) + cI * (TMP2)) + V3[2] * (-cI * (TMP5) + cI * (TMP4)))); V1[3] = denom * (TMP6 * (-cI * (P2[1]) + cI * (P3[1])) + (V2[3] * (-cI * (TMP0) + cI * (TMP2)) + V3[3] * (-cI * (TMP5) + cI * (TMP4)))); V1[4] = denom * (TMP6 * (-cI * (P2[2]) + cI * (P3[2])) + (V2[4] * (-cI * (TMP0) + cI * (TMP2)) + V3[4] * (-cI * (TMP5) + cI * (TMP4)))); V1[5] = denom * (TMP6 * (-cI * (P2[3]) + cI * (P3[3])) + (V2[5] * (-cI * (TMP0) + cI * (TMP2)) + V3[5] * (-cI * (TMP5) + cI * (TMP4)))); } } // end namespace MG5_sm_COLOREA diff --git a/Shower/Dipole/Colorea/read_slha_COLOREA.cc b/Shower/Dipole/Colorea/read_slha_COLOREA.cc --- a/Shower/Dipole/Colorea/read_slha_COLOREA.cc +++ b/Shower/Dipole/Colorea/read_slha_COLOREA.cc @@ -1,144 +1,144 @@ #include #include #include #include "read_slha_COLOREA.h" void SLHABlock_COLOREA::set_entry(std::vector indices, double value) { if (_entries.size() == 0) _indices = indices.size(); else if(indices.size() != _indices) throw "Wrong number of indices in set_entry"; _entries[indices] = value; } double SLHABlock_COLOREA::get_entry(std::vector indices, double def_val) { if (_entries.find(indices) == _entries.end()){ std::cout << "Warning: No such entry in " << _name << ", using default value " << def_val << std::endl; return def_val; } return _entries[indices]; } void SLHAReader_COLOREA::read_slha_file(std::string file_name) { std::ifstream param_card; param_card.open(file_name.c_str(), std::ifstream::in); if(!param_card.good()) throw "Error while opening param card"; std::cout << "\nColorea: Opened slha file " << file_name << " for reading" << std::endl; char buf[200]; std::string line; std::string block(""); while(param_card.good()){ param_card.getline(buf, 200); line = buf; // Change to lowercase transform(line.begin(), line.end(), line.begin(), (int(*)(int)) tolower); if(line != "" && line[0] != '#'){ if(block != ""){ // Look for double index blocks double dindex1, dindex2; double value; std::stringstream linestr2(line); if (linestr2 >> dindex1 >> dindex2 >> value && dindex1 == int(dindex1) and dindex2 == int(dindex2)) { std::vector indices; indices.push_back(int(dindex1)); indices.push_back(int(dindex2)); set_block_entry(block, indices, value); // Done with this line, read next continue; } std::stringstream linestr1(line); // Look for single index blocks if(linestr1 >> dindex1 >> value && dindex1 == int(dindex1)) { std::vector indices; indices.push_back(int(dindex1)); set_block_entry(block, indices, value); // Done with this line, read next continue; } } // Look for block if(line.find("block ") != line.npos){ line = line.substr(6); // Get rid of spaces between block and block name while (line[0] == ' ') line = line.substr(1); // Now find end of block name - int space_pos = line.find(' '); + std::string::size_type space_pos = line.find(' '); if(space_pos != line.npos) line = line.substr(0, space_pos); block = line; continue; } // Look for decay if(line.find("decay ") == 0){ line = line.substr(6); block = ""; std::stringstream linestr(line); int pdg_code; double value; if(linestr >> pdg_code >> value) set_block_entry("decay", pdg_code, value); else std::cout << "Warning: Wrong format for decay block " << line << std::endl; continue; } } } if (_blocks.size() == 0) throw "No information read from SLHA card"; param_card.close(); } double SLHAReader_COLOREA::get_block_entry(std::string block_name, std::vector indices, double def_val) { if (_blocks.find(block_name) == _blocks.end()){ std::cout << "No such block " << block_name << ", using default value " << def_val << std::endl; return def_val; } return _blocks[block_name].get_entry(indices); } double SLHAReader_COLOREA::get_block_entry(std::string block_name, int index, double def_val) { std::vector indices; indices.push_back(index); return get_block_entry(block_name, indices, def_val); } void SLHAReader_COLOREA::set_block_entry(std::string block_name, std::vector indices, double value) { if (_blocks.find(block_name) == _blocks.end()){ SLHABlock_COLOREA block(block_name); _blocks[block_name] = block; } _blocks[block_name].set_entry(indices, value); /* cout << "Set block " << block_name << " entry "; for (int i=0;i < indices.size();i++) cout << indices[i] << " "; cout << "to " << _blocks[block_name].get_entry(indices) << endl;*/ } void SLHAReader_COLOREA::set_block_entry(std::string block_name, int index, double value) { std::vector indices; indices.push_back(index); set_block_entry(block_name, indices, value); } diff --git a/Shower/QTilde/Base/Branching.h b/Shower/QTilde/Base/Branching.h --- a/Shower/QTilde/Base/Branching.h +++ b/Shower/QTilde/Base/Branching.h @@ -1,69 +1,69 @@ // -*- C++ -*- #ifndef HERWIG_Branching_H #define HERWIG_Branching_H // // This is the declaration of the Branching struct. // #include "Herwig/Shower/QTilde/ShowerConfig.h" #include "Herwig/Shower/QTilde/Kinematics/ShowerKinematics.h" namespace Herwig { using namespace ThePEG; /** \ingroup Shower * The branching struct is used to store information on the branching. * The kinematics variable is a pointer to the ShowerKinematics for the branching * The sudakov variable is a pointer to the SudakovFormFactor for the branching * The ids variable is the list of particles in the branching */ struct Branching { /** * Pointer to the ShowerKinematics object for the branching */ ShoKinPtr kinematics; /** * PDG codes of the particles in the branching */ IdList ids; /** * The SudakovFormFactor for the branching */ tSudakovPtr sudakov; /** * The type of radiation line */ ShowerPartnerType type; /** * Whether or not it keep from forced hard emisson */ bool hard; /** * Which of the children is same as incoming */ unsigned int iout; /** * Constructor for the struct * @param a pointer to the ShowerKinematics object for the branching * @param c PDG codes of the particles in the branching * @param d The SudakovFormFactor for the branching */ Branching(ShoKinPtr a, IdList c,tSudakovPtr d,ShowerPartnerType t) : kinematics(a), ids(c), sudakov(d), type(t), hard(false), iout(0) {} /** * Default constructor */ - Branching() : hard(false), iout(0) {} + Branching() : type(ShowerPartnerType::Undefined), hard(false), iout(0) {} }; } #endif /* HERWIG_Branching_H */ diff --git a/Shower/QTilde/Base/ShowerTree.cc b/Shower/QTilde/Base/ShowerTree.cc --- a/Shower/QTilde/Base/ShowerTree.cc +++ b/Shower/QTilde/Base/ShowerTree.cc @@ -1,914 +1,913 @@ // -*- C++ -*- // // ShowerTree.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. // #include "ThePEG/EventRecord/MultiColour.h" #include "ThePEG/Repository/EventGenerator.h" #include "ShowerTree.h" #include "Herwig/Shower/QTilde/Base/ShowerParticle.h" #include "Herwig/Shower/ShowerHandler.h" #include "ThePEG/PDT/DecayMode.h" #include "ThePEG/Handlers/EventHandler.h" #include "ThePEG/Handlers/XComb.h" #include #include "ThePEG/Repository/CurrentGenerator.h" #include "ThePEG/PDT/StandardMatchers.h" using namespace Herwig; using namespace ThePEG; bool ShowerTree::_spaceTime = false; Energy2 ShowerTree::_vmin2 = 0.1*GeV2; namespace { void findBeam(tPPtr & beam, PPtr incoming) { if(beam==incoming) return; while(!beam->children().empty()) { bool found=false; for(unsigned int ix=0;ixchildren().size();++ix) { if(beam->children()[ix]==incoming) { found = true; break; } } if(found) break; beam = beam->children()[0]; } } } ShowerTree::ShowerTree(PerturbativeProcessPtr process) : _parent(), _hasShowered(false) { // get the incoming and outgoing particles and make copies vector original,copy; for(unsigned int ix=0;ixincoming().size();++ix) { original.push_back(process->incoming()[ix].first); copy .push_back(new_ptr(Particle(*original.back()))); } for(unsigned int ix=0;ixoutgoing().size();++ix) { original.push_back(process->outgoing()[ix].first); copy .push_back(new_ptr(Particle(*original.back()))); } // isolate the colour colourIsolate(original,copy); // hard process unsigned int itype(1); if(process->incoming().size()==2) { _wasHard = true; // set the beams and incoming particles tPPair beam = CurrentGenerator::current().currentEvent()->incoming(); findBeam(beam.first ,process->incoming()[0].first); findBeam(beam.second,process->incoming()[1].first); _incoming = make_pair(process->incoming()[0].first, process->incoming()[1].first); double x1(_incoming.first ->momentum().rho()/beam.first ->momentum().rho()); double x2(_incoming.second->momentum().rho()/beam.second->momentum().rho()); // must have two incoming particles assert(_incoming.first && _incoming.second); // set the parent tree _parent=ShowerTreePtr(); for(unsigned int ix=0;ix<2;++ix) { ShowerParticlePtr temp=new_ptr(ShowerParticle(*copy[ix],itype,false)); fixColour(temp); temp->x(ix==0 ? x1 : x2); _incomingLines.insert(make_pair(new_ptr(ShowerProgenitor(original[ix], copy[ix],temp)),temp)); } } // decay process else if(process->incoming().size()==1) { _wasHard = false; itype=2; // create the parent shower particle ShowerParticlePtr sparent(new_ptr(ShowerParticle(*copy[0],itype,false))); fixColour(sparent); _incomingLines.insert(make_pair(new_ptr(ShowerProgenitor(original[0],copy[0],sparent)) ,sparent)); // return if not decayed if(original.size()==1) return; } else assert(false); // create the children assert(copy.size() == original.size()); for (unsigned int ix=process->incoming().size();ixid()==progenitor->id(); } ShowerParticlePtr newpart; if(matches[0]&&matches[1]) { if(parent->showerKinematics()->z()>0.5) newpart=children[0]; else newpart=children[1]; } else if(matches[0]) newpart=children[0]; else if(matches[1]) newpart=children[1]; _outgoingLines[progenitor]=newpart; } void ShowerTree::updateInitialStateShowerProduct(ShowerProgenitorPtr progenitor, ShowerParticlePtr newParent) { _incomingLines[progenitor]=newParent; } void ShowerTree::insertHard(StepPtr pstep, bool ISR, bool) { assert(_incomingLines.size()==2); colourLines().clear(); map::const_iterator cit; // construct the map of colour lines for hard process for(cit=_incomingLines.begin();cit!=_incomingLines.end();++cit) { if(!cit->first->perturbative()) continue; mapColour(cit->first->original(),cit->first->copy()); } map::const_iterator cjt; for(cjt=_outgoingLines.begin();cjt!=_outgoingLines.end();++cjt) { if(!cjt->first->perturbative()) continue; mapColour(cjt->first->original(),cjt->first->copy()); } // initial-state radiation if(ISR) { for(cit=incomingLines().begin();cit!=incomingLines().end();++cit) { ShowerParticlePtr init=(*cit).first->progenitor(); assert(init->thePEGBase()); PPtr original = (*cit).first->original(); if(original->parents().empty()) continue; PPtr hadron= original->parents()[0]; assert(!original->children().empty()); PPtr copy=cit->first->copy(); ParticleVector intermediates=original->children(); for(unsigned int ix=0;ixabandonChild(intermediates[ix]); copy->abandonChild(intermediates[ix]); } // if not from a matrix element correction if(cit->first->perturbative()) { // break mother/daugther relations init->addChild(original); hadron->abandonChild(original); // if particle showers add shower if(cit->first->hasEmitted()) { addInitialStateShower(init,hadron,pstep,false); } // no showering for this particle else { updateColour(init,false); hadron->addChild(init); pstep->addIntermediate(init); init->setLifeLength(Lorentz5Distance()); init->setVertex(LorentzPoint()); } } // from matrix element correction else { // break mother/daugther relations hadron->abandonChild(original); copy->addChild(original); updateColour(copy,false); init->addChild(copy); pstep->addIntermediate(copy); copy->setLifeLength(Lorentz5Distance()); copy->setVertex(LorentzPoint()); // if particle showers add shower if(cit->first->hasEmitted()) { addInitialStateShower(init,hadron,pstep,false); } // no showering for this particle else { updateColour(init,false); hadron->addChild(init); pstep->addIntermediate(init); init->setLifeLength(Lorentz5Distance()); init->setVertex(LorentzPoint()); } } } } else { for(cit=incomingLines().begin();cit!=incomingLines().end();++cit) { ShowerParticlePtr init=(*cit).first->progenitor(); assert(init->thePEGBase()); PPtr original = (*cit).first->original(); if(original->parents().empty()) continue; PPtr hadron= original->parents()[0]; assert(!original->children().empty()); PPtr copy=cit->first->copy(); ParticleVector intermediates=original->children(); for(unsigned int ix=0;ixabandonChild(intermediates[ix]); copy->abandonChild(intermediates[ix]); } // break mother/daugther relations init->addChild(original); hadron->abandonChild(original); // no showering for this particle updateColour(init,false); hadron->addChild(init); pstep->addIntermediate(init); init->setLifeLength(Lorentz5Distance()); init->setVertex(LorentzPoint()); original->setLifeLength(Lorentz5Distance()); original->setVertex(LorentzPoint()); } } // final-state radiation for(cjt=outgoingLines().begin();cjt!=outgoingLines().end();++cjt) { ShowerParticlePtr init=(*cjt).first->progenitor(); assert(init->thePEGBase()); // ZERO the distance of original (*cjt).first->original()->setLifeLength(Lorentz5Distance()); (*cjt).first->original()->setVertex(LorentzPoint()); // if not from a matrix element correction if(cjt->first->perturbative()) { // register the shower particle as a // copy of the one from the hard process tParticleVector parents=init->parents(); for(unsigned int ix=0;ixabandonChild(init); (*cjt).first->original()->addChild(init); pstep->addDecayProduct(init); } // from a matrix element correction else { if(cjt->first->original()==_incoming.first|| cjt->first->original()==_incoming.second) { updateColour((*cjt).first->copy(),false); (*cjt).first->original()->parents()[0]-> addChild((*cjt).first->copy()); pstep->addDecayProduct((*cjt).first->copy()); (*cjt).first->copy()->addChild(init); pstep->addDecayProduct(init); } else { updateColour((*cjt).first->copy(),false); (*cjt).first->original()->addChild((*cjt).first->copy()); pstep->addDecayProduct((*cjt).first->copy()); (*cjt).first->copy()->addChild(init); pstep->addDecayProduct(init); } // ZERO the distance of copy ??? \todo change if add space-time (*cjt).first->copy()->setLifeLength(Lorentz5Distance()); (*cjt).first->copy()->setVertex(LorentzPoint()); } // copy so travels no distance init->setLifeLength(Lorentz5Distance()); init->setVertex(init->parents()[0]->decayVertex()); // sort out the colour updateColour(init,false); // insert shower products addFinalStateShower(init,pstep); } colourLines().clear(); } void ShowerTree::addFinalStateShower(PPtr p, StepPtr s) { // if endpoint assume doesn't travel if(p->children().empty()) { if(p->dataPtr()->stable()||ShowerHandler::currentHandler()->decaysInShower(p->id())) p->setLifeLength(Lorentz5Distance()); else { Energy mass = p->mass()!=ZERO ? p->mass() : p->dataPtr()->mass(); Length ctau = p->dataPtr()->generateLifeTime(mass, p->dataPtr()->width()); Lorentz5Distance lifeLength(ctau,p->momentum().vect()*(ctau/mass)); p->setLifeLength(lifeLength); } return; } // set the space-time distance else { p->setLifeLength(spaceTimeDistance(p)); } ParticleVector::const_iterator child; for(child=p->children().begin(); child != p->children().end(); ++child) { updateColour(*child,false); s->addDecayProduct(*child); (**child).setVertex(p->decayVertex()); addFinalStateShower(*child,s); } } void ShowerTree::addInitialStateShower(PPtr p, PPtr hadron, StepPtr s, bool addchildren) { // Each parton here should only have one parent if(!p->parents().empty()) { if(p->parents().size()!=1) throw Exception() << "Particle must only have one parent in ShowerTree" << "::addInitialStateShower" << Exception::runerror; // set the space-time distances if(addchildren) { p->setLifeLength(spaceTimeDistance(p)); p->setVertex(p->children()[0]->vertex()-p->lifeLength()); } else { p->setLifeLength(spaceTimeDistance(p)); p->setVertex(-p->lifeLength()); } // recurse addInitialStateShower(p->parents()[0],hadron,s); } else { hadron->addChild(p); s->addIntermediate(p); p->setVertex(p->children()[0]->vertex()); p->setLifeLength(Lorentz5Distance()); } updateColour(p,false); // if not adding children return if(!addchildren) return; // add children ParticleVector::const_iterator child; for(child = p->children().begin(); child != p->children().end(); ++child) { // if a final-state particle update the colour ShowerParticlePtr schild = dynamic_ptr_cast(*child); (**child).setVertex(p->decayVertex()); if(schild && schild->isFinalState()) updateColour(*child,false); // if there are grandchildren of p if(!(*child)->children().empty()) { // Add child as intermediate s->addIntermediate(*child); // If child is shower particle and final-state, add children if(schild && schild->isFinalState()) addFinalStateShower(schild,s); } else s->addDecayProduct(*child); } } void ShowerTree::update(PerturbativeProcessPtr newProcess) { // must be one incoming particle assert(_incomingLines.size()==1); colourLines().clear(); // copy the particles and isolate the colour vector original,copy; for(unsigned int ix=0;ixincoming().size();++ix) { original.push_back(newProcess->incoming()[ix].first); copy .push_back(new_ptr(Particle(*original.back()))); } for(unsigned int ix=0;ixoutgoing().size();++ix) { original.push_back(newProcess->outgoing()[ix].first); copy .push_back(new_ptr(Particle(*original.back()))); } // isolate the colour colourIsolate(original,copy); // make the new progenitor ShowerParticlePtr stemp=new_ptr(ShowerParticle(*copy[0],2,false)); fixColour(stemp); ShowerProgenitorPtr newprog=new_ptr(ShowerProgenitor(original[0],copy[0],stemp)); _incomingLines.clear(); _incomingLines.insert(make_pair(newprog,stemp)); // create the children assert(copy.size() == original.size()); for (unsigned int ix=newProcess->incoming().size();ix_treelinks.empty()) final = _parent->_treelinks[this].second; else final=_incomingLines.begin()->first->original(); // construct the map of colour lines PPtr copy=_incomingLines.begin()->first->copy(); mapColour(final,copy); // now this is the ONE instance of the particle which should have a life length // \todo change if space-time picture added // set the lifelength, need this so that still in right direction after // any ISR recoils Length ctau = copy->lifeTime(); Lorentz5Distance lifeLength(ctau,final->momentum().vect()*(ctau/final->mass())); final->setLifeLength(lifeLength); // initial-state radiation if(ISR&&!_incomingLines.begin()->first->progenitor()->children().empty()) { ShowerParticlePtr init=_incomingLines.begin()->first->progenitor(); updateColour(init,false); final->addChild(init); pstep->addDecayProduct(init); // just a copy doesn't travel init->setLifeLength(Lorentz5Distance()); init->setVertex(final->decayVertex()); // insert shower products addFinalStateShower(init,pstep); // sort out colour final=_incomingLines.begin()->second; colourLines().clear(); mapColour(final,copy); } // get the decaying particles // make the copy tColinePair cline=make_pair(copy->colourLine(),copy->antiColourLine()); updateColour(copy,false); // sort out sinks and sources if needed if(cline.first) { if(cline.first->sourceNeighbours().first) { copy->colourLine()->setSourceNeighbours(cline.first->sourceNeighbours().first, cline.first->sourceNeighbours().second); } else if (cline.first->sinkNeighbours().first) { copy->colourLine()->setSinkNeighbours(cline.first->sinkNeighbours().first, cline.first->sinkNeighbours().second); } } if(cline.second) { if(cline.second->sourceNeighbours().first) { copy->antiColourLine()->setSourceNeighbours(cline.second->sourceNeighbours().first, cline.second->sourceNeighbours().second); } else if (cline.second->sinkNeighbours().first) { copy->antiColourLine()->setSinkNeighbours(cline.second->sinkNeighbours().first, cline.second->sinkNeighbours().second); } } // copy of the one from the hard process tParticleVector dpar=copy->parents(); for(unsigned int ix=0;ixabandonChild(copy); final->addChild(copy); pstep->addDecayProduct(copy); // just a copy does not move copy->setLifeLength(Lorentz5Distance()); copy->setVertex(final->decayVertex()); // final-state radiation map::const_iterator cit; for(cit=outgoingLines().begin();cit!=outgoingLines().end();++cit) { ShowerParticlePtr init=cit->first->progenitor(); // ZERO the distance init->setLifeLength(Lorentz5Distance()); if(!init->thePEGBase()) throw Exception() << "Final-state particle must have a ThePEGBase" << " in ShowerTree::insertDecay()" << Exception::runerror; // if not from matrix element correction if(cit->first->perturbative()) { // add the child updateColour(cit->first->copy(),false); PPtr orig=cit->first->original(); orig->setLifeLength(Lorentz5Distance()); orig->setVertex(copy->decayVertex()); copy->addChild(orig); pstep->addDecayProduct(orig); orig->addChild(cit->first->copy()); pstep->addDecayProduct(cit->first->copy()); // register the shower particle as a // copy of the one from the hard process tParticleVector parents=init->parents(); for(unsigned int ix=0;ixabandonChild(init);} (*cit).first->copy()->addChild(init); pstep->addDecayProduct(init); updateColour(init,false); } // from a matrix element correction else { if(copy->children().end()== find(copy->children().begin(),copy->children().end(), cit->first->original())) { updateColour(cit->first->original(),false); copy->addChild(cit->first->original()); pstep->addDecayProduct(cit->first->original()); } updateColour(cit->first->copy(),false); cit->first->original()->addChild(cit->first->copy()); pstep->addDecayProduct(cit->first->copy()); // register the shower particle as a // copy of the one from the hard process tParticleVector parents=init->parents(); for(unsigned int ix=0;ixabandonChild(init);} (*cit).first->copy()->addChild(init); pstep->addDecayProduct(init); updateColour(init,false); } // ZERO the distances as just copies cit->first->copy()->setLifeLength(Lorentz5Distance()); init->setLifeLength(Lorentz5Distance()); cit->first->copy()->setVertex(copy->decayVertex()); init->setVertex(copy->decayVertex()); // insert shower products addFinalStateShower(init,pstep); } colourLines().clear(); } void ShowerTree::clear() { // reset the has showered flag _hasShowered=false; // clear the colour map colourLines().clear(); map::const_iterator cit; map::const_iterator cjt; // abandon the children of the outgoing particles for(cit=_outgoingLines.begin();cit!=_outgoingLines.end();++cit) { ShowerParticlePtr orig=cit->first->progenitor(); orig->set5Momentum(cit->first->copy()->momentum()); ParticleVector children=orig->children(); for(unsigned int ix=0;ixabandonChild(children[ix]); _outgoingLines[cit->first]=orig; cit->first->hasEmitted(false); cit->first->reconstructed(ShowerProgenitor::notReconstructed); } // forward products _forward.clear(); for(cit=_outgoingLines.begin();cit!=_outgoingLines.end();++cit) _forward.insert(cit->first->progenitor()); // if a decay if(!_wasHard) { ShowerParticlePtr orig=_incomingLines.begin()->first->progenitor(); orig->set5Momentum(_incomingLines.begin()->first->copy()->momentum()); ParticleVector children=orig->children(); for(unsigned int ix=0;ixabandonChild(children[ix]); _incomingLines.begin()->first->reconstructed(ShowerProgenitor::notReconstructed); } // if a hard process else { for(cjt=_incomingLines.begin();cjt!=_incomingLines.end();++cjt) { tPPtr parent = cjt->first->original()->parents().empty() ? tPPtr() : cjt->first->original()->parents()[0]; ShowerParticlePtr temp= new_ptr(ShowerParticle(*cjt->first->copy(), cjt->first->progenitor()->perturbative(), cjt->first->progenitor()->isFinalState())); fixColour(temp); temp->x(cjt->first->progenitor()->x()); cjt->first->hasEmitted(false); if(!(cjt->first->progenitor()==cjt->second)&&cjt->second&&parent) parent->abandonChild(cjt->second); cjt->first->progenitor(temp); _incomingLines[cjt->first]=temp; cjt->first->reconstructed(ShowerProgenitor::notReconstructed); } } clearTransforms(); } void ShowerTree::resetShowerProducts() { map::const_iterator cit; map::const_iterator cjt; _forward.clear(); for(cjt=_outgoingLines.begin();cjt!=_outgoingLines.end();++cjt) _forward.insert(cjt->second); } void ShowerTree::updateAfterShower(ShowerDecayMap & decay) { // update the links map::const_iterator mit; map >::iterator tit; for(tit=_treelinks.begin();tit!=_treelinks.end();++tit) { if(tit->second.first) { mit=_outgoingLines.find(tit->second.first); if(mit!=_outgoingLines.end()) tit->second.second=mit->second; } } // get the particles coming from those in the hard process set hard; for(mit=_outgoingLines.begin();mit!=_outgoingLines.end();++mit) hard.insert(mit->second); // find the shower particles which should be decayed in the // shower but didn't come from the hard process set::const_iterator cit; // TODO construct _forward here? for(cit=_forward.begin();cit!=_forward.end();++cit) { if((ShowerHandler::currentHandler()->decaysInShower((**cit).id())&& !(**cit).dataPtr()->stable()) && hard.find(*cit)==hard.end()) { PerturbativeProcessPtr newProcess(new_ptr(PerturbativeProcess())); newProcess->incoming().push_back(make_pair(*cit,PerturbativeProcessPtr())); ShowerTreePtr newtree=new_ptr(ShowerTree(newProcess)); newtree->setParents(); newtree->_parent=this; Energy width=(**cit).dataPtr()->generateWidth((**cit).mass()); decay.insert(make_pair(width,newtree)); _treelinks.insert(make_pair(newtree, make_pair(tShowerProgenitorPtr(),*cit))); } } } void ShowerTree::addFinalStateBranching(ShowerParticlePtr parent, const ShowerParticleVector & children) { assert(children.size()==2); _forward.erase(parent); for(unsigned int ix=0; ix >::const_iterator tit; for(tit=_treelinks.begin();tit!=_treelinks.end();++tit) tit->first->_parent=this; } vector ShowerTree::extractProgenitors() { // extract the particles from the ShowerTree map::const_iterator mit; vector ShowerHardJets; for(mit=incomingLines().begin();mit!=incomingLines().end();++mit) ShowerHardJets.push_back((*mit).first); map::const_iterator mjt; for(mjt=outgoingLines().begin();mjt!=outgoingLines().end();++mjt) ShowerHardJets.push_back((*mjt).first); return ShowerHardJets; } void ShowerTree::transform(const LorentzRotation & boost, bool applyNow) { if(applyNow) { // now boost all the particles map::const_iterator cit; // incoming for(cit=_incomingLines.begin();cit!=_incomingLines.end();++cit) { cit->first->progenitor()->deepTransform(boost); cit->first->copy()->deepTransform(boost); } // outgoing map::const_iterator cjt; for(cjt=_outgoingLines.begin();cjt!=_outgoingLines.end();++cjt) { cjt->first->progenitor()->deepTransform(boost); cjt->first->copy()->deepTransform(boost); } } else { _transforms.transform(boost); } // child trees for(map >::const_iterator tit=_treelinks.begin();tit!=_treelinks.end();++tit) tit->first->transform(boost,applyNow); } void ShowerTree::applyTransforms() { // now boost all the particles map::const_iterator cit; // incoming for(cit=_incomingLines.begin();cit!=_incomingLines.end();++cit) { cit->first->progenitor()->deepTransform(_transforms); cit->first->copy()->deepTransform(_transforms); } // outgoing map::const_iterator cjt; for(cjt=_outgoingLines.begin();cjt!=_outgoingLines.end();++cjt) { cjt->first->progenitor()->deepTransform(_transforms); cjt->first->copy()->deepTransform(_transforms); } // child trees for(map >::const_iterator tit=_treelinks.begin();tit!=_treelinks.end();++tit) tit->first->applyTransforms(); _transforms = LorentzRotation(); } void ShowerTree::clearTransforms() { _transforms = LorentzRotation(); // // child trees // for(map >::const_iterator // tit=_treelinks.begin();tit!=_treelinks.end();++tit) // tit->first->clearTransforms(); } void ShowerTree::fixColour(tShowerParticlePtr part) { if(!part->colourInfo()->colourLines().empty()) { if(part->colourInfo()->colourLines().size()==1) { ColinePtr line=part->colourLine(); if(line) { line->removeColoured(part); line->addColoured(part); } } else { Ptr::pointer colour = dynamic_ptr_cast::pointer>(part->colourInfo()); vector lines = colour->colourLines(); for(unsigned int ix=0;ix(lines[ix]); if(line) { line->removeColoured(part); line->addColoured(part); } } } } if(!part->colourInfo()->antiColourLines().empty()) { if(part->colourInfo()->antiColourLines().size()==1) { ColinePtr line=part->antiColourLine(); if(line) { line->removeAntiColoured(part); line->addAntiColoured(part); } } else { Ptr::pointer colour = dynamic_ptr_cast::pointer>(part->colourInfo()); vector lines = colour->antiColourLines(); for(unsigned int ix=0;ix(lines[ix]); if(line) { line->removeAntiColoured(part); line->addAntiColoured(part); } } } } } vector ShowerTree::extractProgenitorParticles() { vector particles; // incoming particles for(map::const_iterator cit=incomingLines().begin(); cit!=incomingLines().end();++cit) particles.push_back(cit->first->progenitor()); // outgoing particles for(map::const_iterator cjt=outgoingLines().begin(); cjt!=outgoingLines().end();++cjt) particles.push_back(cjt->first->progenitor()); return particles; } Lorentz5Distance ShowerTree::spaceTimeDistance(tPPtr particle) { if(!_spaceTime) return Lorentz5Distance(); Energy2 q2 = particle->mass() > ZERO ? sqr(particle->mass()) : -sqr(particle->mass()); const tcPDPtr data = particle->dataPtr(); // calculate width imposing min value Energy conMass = max(data->constituentMass(),200*MeV); Energy width = max(data->generateWidth(particle->mass()),_vmin2/conMass); // offshellness Energy2 offShell = q2-sqr(data->constituentMass()); if(abs(offShell)<1e-10*GeV2) offShell = ZERO; InvEnergy2 fact = UseRandom::rndExp(1./sqrt((sqr(offShell)+sqr(q2*width/conMass)))); Lorentz5Distance output = (hbarc*fact)*particle->momentum(); return output; } void ShowerTree::constructTrees(ShowerTreePtr & hardTree, ShowerDecayMap & decayTrees, PerturbativeProcessPtr hard, DecayProcessMap decay) { map treeMap; // convert the hard process if(hardTree) { if(hardTree->isDecay()) hardTree->update(hard); } else { hardTree = new_ptr(ShowerTree(hard)); } treeMap.insert(make_pair(hard,hardTree)); for(DecayProcessMap::const_iterator it=decay.begin();it!=decay.end();++it) { ShowerTreePtr newTree = new_ptr(ShowerTree(it->second)); treeMap.insert(make_pair(it->second,newTree)); decayTrees.insert(make_pair(it->first,newTree)); } // set up the links between the trees for(map::const_iterator it=treeMap.begin(); it!=treeMap.end();++it) { // links to daughter trees map::const_iterator mit; for(mit=it->second->_outgoingLines.begin(); mit!=it->second->_outgoingLines.end();++mit) { unsigned int iloc=0; for(;ilocfirst->outgoing().size();++iloc) { if(it->first->outgoing()[iloc].first==mit->first->original()) break; } if(it->first->outgoing()[iloc].second) it->second->_treelinks.insert(make_pair(treeMap[it->first->outgoing()[iloc].second], make_pair(mit->first,mit->first->progenitor()))); } // links to parent trees if(it->first->incoming()[0].second) { it->second->_parent = treeMap[it->first->incoming()[0].second]; } } } namespace { Lorentz5Momentum sumMomenta(tPPtr particle) { if(particle->children().empty()) return particle->momentum(); Lorentz5Momentum output; for(unsigned int ix=0;ixchildren().size();++ix) { output += sumMomenta(particle->children()[ix]); } return output; } } void ShowerTree::checkMomenta() { vector pin; for(map::const_iterator it=_incomingLines.begin(); it!=_incomingLines.end();++it) { if(isDecay()) { tPPtr parent = it->first->progenitor(); pin.push_back(parent->momentum()); while(!parent->children().empty()) { pin.back() -= sumMomenta(parent->children()[1]); parent = parent->children()[0]; } } else { tPPtr parent = it->second; pin.push_back(parent->momentum()); while(!parent->children().empty()&&parent->children().size()==2) { pin.back() -= sumMomenta(parent->children()[1]); parent = parent->children()[0]; if(parent==it->first->progenitor()) break; } } } vector pout; for(map::const_iterator it= _outgoingLines.begin(); it!=_outgoingLines.end();++it) { pout.push_back(sumMomenta(it->first->progenitor())); } Lorentz5Momentum psum; for(unsigned int ix=0;ix< pin.size();++ix) { CurrentGenerator::log() << "pin " << ix << pin[ix]/GeV << "\n"; psum +=pin[ix]; } CurrentGenerator::log() << "In total " << psum/GeV << " " << psum.m()/GeV << "\n"; Lorentz5Momentum psum2; for(unsigned int ix=0;ix< pout.size();++ix) { CurrentGenerator::log() << "pout " << ix << pout[ix]/GeV << "\n"; psum2 +=pout[ix]; } CurrentGenerator::log() << "Out total " << psum2/GeV << " " << psum2.m()/GeV << "\n"; CurrentGenerator::log() << "Total " << (psum-psum2)/GeV << "\n"; } RealEmissionProcessPtr ShowerTree::perturbativeProcess() { RealEmissionProcessPtr output(new_ptr(RealEmissionProcess())); // add the incoming particles unsigned int ix=0; pair x; for(map::const_iterator it=_incomingLines.begin(); it!=_incomingLines.end();++it) { output->bornIncoming().push_back(it->first->progenitor()); if(!it->first->original()->parents().empty()) output->hadrons().push_back(it->first->original()->parents()[0]); else output->hadrons().push_back(it->first->progenitor()); if(ix==0) x.first = it->second->x(); else x.second = it->second->x(); ++ix; } output->x(x); // add the outgoing particles for(map::const_iterator it= _outgoingLines.begin(); it!=_outgoingLines.end();++it) { output->bornOutgoing().push_back(it->first->progenitor()); } return output; } void ShowerTree::setVetoes(const map & pTs, unsigned int type) { if(type==1||type==3) { for(map::const_iterator it=_incomingLines.begin(); it!=_incomingLines.end();++it) { for(map::const_iterator jt=pTs.begin();jt!=pTs.end();++jt) it->first->maximumpT(jt->second,jt->first); } } if(type==2||type==3) { for(map::const_iterator it= _outgoingLines.begin(); it!=_outgoingLines.end();++it) { for(map::const_iterator jt=pTs.begin();jt!=pTs.end();++jt) it->first->maximumpT(jt->second,jt->first); } } } diff --git a/Shower/QTilde/Kinematics/KinematicsReconstructor.cc b/Shower/QTilde/Kinematics/KinematicsReconstructor.cc --- a/Shower/QTilde/Kinematics/KinematicsReconstructor.cc +++ b/Shower/QTilde/Kinematics/KinematicsReconstructor.cc @@ -1,2944 +1,2945 @@ // -*- C++ -*- // // KinematicsReconstructor.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 KinematicsReconstructor class. // #include "KinematicsReconstructor.h" #include "ThePEG/PDT/EnumParticles.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/EventRecord/Event.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/RefVector.h" #include "Herwig/Shower/QTilde/Base/PartnerFinder.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/Shower/QTilde/SplittingFunctions/SplittingFunction.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/EventRecord/ColourLine.h" #include "ThePEG/Utilities/DescribeClass.h" #include "Herwig/Shower/QTilde/QTildeShowerHandler.h" #include +#include "KinematicsReconstructor.tcc" using namespace Herwig; DescribeClass describeKinematicsReconstructor("Herwig::KinematicsReconstructor", "HwShower.so"); namespace { /** * Struct to order the jets in off-shellness */ struct JetOrdering { bool operator() (const JetKinStruct & j1, const JetKinStruct & j2) { Energy diff1 = j1.q.m()-j1.p.m(); Energy diff2 = j2.q.m()-j2.p.m(); if(diff1!=diff2) { return diff1>diff2; } else if( j1.q.e() != j2.q.e() ) return j1.q.e()>j2.q.e(); else return j1.parent->uniqueId>j2.parent->uniqueId; } }; } void KinematicsReconstructor::persistentOutput(PersistentOStream & os) const { os << _reconopt << _initialBoost << ounit(_minQ,GeV) << _noRescale << _noRescaleVector << _finalStateReconOption << _initialStateReconOption; } void KinematicsReconstructor::persistentInput(PersistentIStream & is, int) { is >> _reconopt >> _initialBoost >> iunit(_minQ,GeV) >> _noRescale >> _noRescaleVector >> _finalStateReconOption >> _initialStateReconOption; } void KinematicsReconstructor::Init() { static ClassDocumentation documentation ( "This class is responsible for the kinematics reconstruction of the showering,", " including the kinematics reshuffling necessary to compensate for the recoil" "of the emissions." ); static Switch interfaceReconstructionOption ("ReconstructionOption", "Option for the kinematics reconstruction", &KinematicsReconstructor::_reconopt, 0, false, false); static SwitchOption interfaceReconstructionOptionGeneral (interfaceReconstructionOption, "General", "Use the general solution which ignores the colour structure for all processes", 0); static SwitchOption interfaceReconstructionOptionColour (interfaceReconstructionOption, "Colour", "Use the colour structure of the process to determine the reconstruction procedure.", 1); static SwitchOption interfaceReconstructionOptionColour2 (interfaceReconstructionOption, "Colour2", "Make the most use possible of the colour structure of the process to determine the reconstruction procedure. " "Start with FF, then IF then II colour connections", 2); static SwitchOption interfaceReconstructionOptionColour3 (interfaceReconstructionOption, "Colour3", "Make the most use possible of the colour structure of the process to determine the reconstruction procedure. " "Do the colour connections in order of the pT's emitted in the shower starting with the hardest." " The colour partner is fully reconstructed at the same time.", 3); static SwitchOption interfaceReconstructionOptionColour4 (interfaceReconstructionOption, "Colour4", "Make the most use possible of the colour structure of the process to determine the reconstruction procedure. " "Do the colour connections in order of the pT's emitted in the shower starting with the hardest, while leaving" " the colour partner on mass-shell", 4); static Parameter interfaceMinimumQ2 ("MinimumQ2", "The minimum Q2 for the reconstruction of initial-final systems", &KinematicsReconstructor::_minQ, GeV, 0.001*GeV, 1e-6*GeV, 10.0*GeV, false, false, Interface::limited); static RefVector interfaceNoRescale ("NoRescale", "Particles which shouldn't be rescaled to be on shell by the shower", &KinematicsReconstructor::_noRescaleVector, -1, false, false, true, false, false); static Switch interfaceInitialInitialBoostOption ("InitialInitialBoostOption", "Option for how the boost from the system before ISR to that after ISR is applied.", &KinematicsReconstructor::_initialBoost, 0, false, false); static SwitchOption interfaceInitialInitialBoostOptionOneBoost (interfaceInitialInitialBoostOption, "OneBoost", "Apply one boost from old CMS to new CMS", 0); static SwitchOption interfaceInitialInitialBoostOptionLongTransBoost (interfaceInitialInitialBoostOption, "LongTransBoost", "First apply a longitudinal and then a transverse boost", 1); static Switch interfaceFinalStateReconOption ("FinalStateReconOption", "Option for how to reconstruct the momenta of the final-state system", &KinematicsReconstructor::_finalStateReconOption, 0, false, false); static SwitchOption interfaceFinalStateReconOptionDefault (interfaceFinalStateReconOption, "Default", "All the momenta are rescaled in the rest frame", 0); static SwitchOption interfaceFinalStateReconOptionMostOffShell (interfaceFinalStateReconOption, "MostOffShell", "All particles put on the new-mass shell and then the most off-shell and" " recoiling system are rescaled to ensure 4-momentum is conserved.", 1); static SwitchOption interfaceFinalStateReconOptionRecursive (interfaceFinalStateReconOption, "Recursive", "Recursively put on shell by putting the most off-shell particle which" " hasn't been rescaled on-shell by rescaling the particles and the recoiling system. ", 2); static SwitchOption interfaceFinalStateReconOptionRestMostOffShell (interfaceFinalStateReconOption, "RestMostOffShell", "The most off-shell is put on shell by rescaling it and the recoiling system," " the recoiling system is then put on-shell in its rest frame.", 3); static SwitchOption interfaceFinalStateReconOptionRestRecursive (interfaceFinalStateReconOption, "RestRecursive", "As 3 but recursive treated the currently most-off shell," " only makes a difference if more than 3 partons.", 4); static Switch interfaceInitialStateReconOption ("InitialStateReconOption", "Option for the reconstruction of initial state radiation", &KinematicsReconstructor::_initialStateReconOption, 0, false, false); static SwitchOption interfaceInitialStateReconOptionRapidity (interfaceInitialStateReconOption, "Rapidity", "Preserve shat and rapidity", 0); static SwitchOption interfaceInitialStateReconOptionLongitudinal (interfaceInitialStateReconOption, "Longitudinal", "Preserve longitudinal momentum", 1); static SwitchOption interfaceInitialStateReconOptionSofterFraction (interfaceInitialStateReconOption, "SofterFraction", "Preserve the momentum fraction of the parton which has emitted softer.", 2); } void KinematicsReconstructor::doinit() { Interfaced::doinit(); _noRescale = set(_noRescaleVector.begin(),_noRescaleVector.end()); } bool KinematicsReconstructor:: reconstructTimeLikeJet(const tShowerParticlePtr particleJetParent) const { assert(particleJetParent); bool emitted=true; // if this is not a fixed point in the reconstruction if( !particleJetParent->children().empty() ) { // if not a reconstruction fixpoint, dig deeper for all children: for ( ParticleVector::const_iterator cit = particleJetParent->children().begin(); cit != particleJetParent->children().end(); ++cit ) reconstructTimeLikeJet(dynamic_ptr_cast(*cit)); } // it is a reconstruction fixpoint, ie kinematical data has to be available else { // check if the parent was part of the shower ShowerParticlePtr jetGrandParent; if(!particleJetParent->parents().empty()) jetGrandParent= dynamic_ptr_cast (particleJetParent->parents()[0]); // update if so if (jetGrandParent) { if (jetGrandParent->showerKinematics()) { if(particleJetParent->id()==_progenitor->id()&& !_progenitor->data().stable()&&abs(_progenitor->data().id())!=ParticleID::tauminus) { jetGrandParent->showerKinematics()->reconstructLast(particleJetParent, _progenitor->mass()); } else { jetGrandParent->showerKinematics()->reconstructLast(particleJetParent); } } } // otherwise else { Energy dm = ShowerHandler::currentHandler()->retConstituentMasses()? particleJetParent->data().constituentMass(): particleJetParent->data().mass(); if (abs(dm-particleJetParent->momentum().m())>0.001*MeV &&(particleJetParent->dataPtr()->stable() || abs(particleJetParent->id())==ParticleID::tauminus) &&particleJetParent->id()!=ParticleID::gamma &&_noRescale.find(particleJetParent->dataPtr())==_noRescale.end()) { Lorentz5Momentum dum = particleJetParent->momentum(); dum.setMass(dm); dum.rescaleEnergy(); particleJetParent->set5Momentum(dum); } else { emitted=false; } } } // recursion has reached an endpoint once, ie we can reconstruct the // kinematics from the children. if( !particleJetParent->children().empty() ) particleJetParent->showerKinematics() ->reconstructParent( particleJetParent, particleJetParent->children() ); return emitted; } bool KinematicsReconstructor:: reconstructHardJets(ShowerTreePtr hard, const map > & intrinsic, ShowerInteraction type, bool switchRecon) const { _currentTree = hard; _intrinsic=intrinsic; // extract the particles from the ShowerTree vector ShowerHardJets=hard->extractProgenitors(); for(unsigned int ix=0;ixprogenitor()] = vector(); } for(map >::const_iterator tit = _currentTree->treelinks().begin(); tit != _currentTree->treelinks().end();++tit) { _treeBoosts[tit->first] = vector(); } try { // old recon method, using new member functions if(_reconopt == 0 || switchRecon ) { reconstructGeneralSystem(ShowerHardJets); } // reconstruction based on coloured systems else if( _reconopt == 1) { reconstructColourSinglets(ShowerHardJets,type); } // reconstruction of FF, then IF, then II else if( _reconopt == 2) { reconstructFinalFirst(ShowerHardJets); } // reconstruction based on coloured systems else if( _reconopt == 3 || _reconopt == 4) { reconstructColourPartner(ShowerHardJets); } else assert(false); } catch(KinematicsReconstructionVeto) { _progenitor=tShowerParticlePtr(); _intrinsic.clear(); for(map >::const_iterator bit=_boosts.begin();bit!=_boosts.end();++bit) { for(vector::const_reverse_iterator rit=bit->second.rbegin();rit!=bit->second.rend();++rit) { LorentzRotation rot = rit->inverse(); bit->first->transform(rot); } } _boosts.clear(); for(map >::const_iterator bit=_treeBoosts.begin();bit!=_treeBoosts.end();++bit) { for(vector::const_reverse_iterator rit=bit->second.rbegin();rit!=bit->second.rend();++rit) { LorentzRotation rot = rit->inverse(); bit->first->transform(rot,false); } } _currentTree = tShowerTreePtr(); _treeBoosts.clear(); return false; } catch (Exception & ex) { _progenitor=tShowerParticlePtr(); _intrinsic.clear(); _currentTree = tShowerTreePtr(); _boosts.clear(); _treeBoosts.clear(); throw ex; } _progenitor=tShowerParticlePtr(); _intrinsic.clear(); // ensure x<1 for(map::const_iterator cit=hard->incomingLines().begin();cit!=hard->incomingLines().end();++cit) { tPPtr parent = cit->first->progenitor(); while (!parent->parents().empty()) { parent = parent->parents()[0]; } tPPtr hadron; if ( cit->first->original()->parents().empty() ) { hadron = cit->first->original(); } else { hadron = cit->first->original()->parents()[0]; } if( ! (hadron->id() == parent->id() && hadron->children().size() <= 1) && parent->momentum().rho() > hadron->momentum().rho()) { _progenitor=tShowerParticlePtr(); _intrinsic.clear(); for(map >::const_iterator bit=_boosts.begin();bit!=_boosts.end();++bit) { for(vector::const_reverse_iterator rit=bit->second.rbegin();rit!=bit->second.rend();++rit) { LorentzRotation rot = rit->inverse(); bit->first->transform(rot); } } _boosts.clear(); for(map >::const_iterator bit=_treeBoosts.begin();bit!=_treeBoosts.end();++bit) { for(vector::const_reverse_iterator rit=bit->second.rbegin();rit!=bit->second.rend();++rit) { LorentzRotation rot = rit->inverse(); bit->first->transform(rot,false); } } _currentTree = tShowerTreePtr(); _treeBoosts.clear(); return false; } } _boosts.clear(); _treeBoosts.clear(); _currentTree = tShowerTreePtr(); return true; } double KinematicsReconstructor::solveKfactor(const Energy & root_s, const JetKinVect & jets) const { Energy2 s = sqr(root_s); // must be at least two jets if ( jets.size() < 2) throw KinematicsReconstructionVeto(); // sum of jet masses must be less than roots if(momConsEq( 0.0, root_s, jets )>ZERO) throw KinematicsReconstructionVeto(); // if two jets simple solution if ( jets.size() == 2 ) { static const Energy2 eps = 1.0e-4 * MeV2; if ( sqr(jets[0].p.x()+jets[1].p.x()) < eps && sqr(jets[0].p.y()+jets[1].p.y()) < eps && sqr(jets[0].p.z()+jets[1].p.z()) < eps ) { Energy test = (jets[0].p+jets[1].p).vect().mag(); if(test > 1.0e-4 * MeV) throw KinematicsReconstructionVeto(); if ( jets[0].p.vect().mag2() < eps ) throw KinematicsReconstructionVeto(); Energy2 m1sq(jets[0].q.m2()),m2sq(jets[1].q.m2()); return sqrt( ( sqr(s - m1sq - m2sq) - 4.*m1sq*m2sq ) /(4.*s*jets[0].p.vect().mag2()) ); } else throw KinematicsReconstructionVeto(); } // i.e. jets.size() > 2, numerically // check convergence, if it's a problem maybe use Newton iteration? else { double k1 = 0.,k2 = 1.,k = 0.; if ( momConsEq( k1, root_s, jets ) < ZERO ) { while ( momConsEq( k2, root_s, jets ) < ZERO ) { k1 = k2; k2 *= 2; } while ( fabs( (k1 - k2)/(k1 + k2) ) > 1.e-10 ) { if( momConsEq( k2, root_s, jets ) == ZERO ) { return k2; } else { k = (k1+k2)/2.; if ( momConsEq( k, root_s, jets ) > ZERO ) { k2 = k; } else { k1 = k; } } } return k1; } else throw KinematicsReconstructionVeto(); } throw KinematicsReconstructionVeto(); } bool KinematicsReconstructor:: reconstructSpaceLikeJet( const tShowerParticlePtr p) const { bool emitted = true; tShowerParticlePtr child; tShowerParticlePtr parent; if(!p->parents().empty()) parent = dynamic_ptr_cast(p->parents()[0]); if(parent) { emitted=true; reconstructSpaceLikeJet(parent); } // if branching reconstruct time-like child if(p->children().size()==2) child = dynamic_ptr_cast(p->children()[1]); if(p->perturbative()==0 && child) { dynamic_ptr_cast(p->children()[0])-> showerKinematics()->reconstructParent(p,p->children()); if(!child->children().empty()) { _progenitor=child; reconstructTimeLikeJet(child); // calculate the momentum of the particle Lorentz5Momentum pnew=p->momentum()-child->momentum(); pnew.rescaleMass(); p->children()[0]->set5Momentum(pnew); } } return emitted; } Boost KinematicsReconstructor:: solveBoostBeta( const double k, const Lorentz5Momentum & newq, const Lorentz5Momentum & oldp ) { // try something different, purely numerical first: // a) boost to rest frame of newq, b) boost with kp/E Energy q = newq.vect().mag(); Energy2 qs = sqr(q); Energy2 Q2 = newq.m2(); Energy kp = k*(oldp.vect().mag()); Energy2 kps = sqr(kp); // usually we take the minus sign, since this boost will be smaller. // we only require |k \vec p| = |\vec q'| which leaves the sign of // the boost open but the 'minus' solution gives a smaller boost // parameter, i.e. the result should be closest to the previous // result. this is to be changed if we would get many momentum // conservation violations at the end of the shower from a hard // process. double betam = (q*sqrt(qs + Q2) - kp*sqrt(kps + Q2))/(kps + qs + Q2); // move directly to 'return' Boost beta = -betam*(k/kp)*oldp.vect(); // note that (k/kp)*oldp.vect() = oldp.vect()/oldp.vect().mag() but cheaper. // leave this out if it's running properly! if ( betam >= 0 ) return beta; else return Boost(0., 0., 0.); } bool KinematicsReconstructor:: reconstructDecayJets(ShowerTreePtr decay, ShowerInteraction) const { _currentTree = decay; // extract the particles from the ShowerTree vector ShowerHardJets=decay->extractProgenitors(); for(unsigned int ix=0;ixprogenitor()] = vector(); } for(map >::const_iterator tit = _currentTree->treelinks().begin(); tit != _currentTree->treelinks().end();++tit) { _treeBoosts[tit->first] = vector(); } try { bool radiated[2]={false,false}; // find the decaying particle and check if particles radiated ShowerProgenitorPtr initial; for(unsigned int ix=0;ixprogenitor()->isFinalState()) { radiated[1] |=ShowerHardJets[ix]->hasEmitted(); } else { initial=ShowerHardJets[ix]; radiated[0]|=ShowerHardJets[ix]->hasEmitted(); } } // find boost to the rest frame if needed Boost boosttorest=-initial->progenitor()->momentum().boostVector(); double gammarest = initial->progenitor()->momentum().e()/ initial->progenitor()->momentum().mass(); // check if need to boost to rest frame bool gottaBoost = (boosttorest.mag() > 1e-12); // if initial state radiation reconstruct the jet and set up the basis vectors Lorentz5Momentum pjet; Lorentz5Momentum nvect; // find the partner ShowerParticlePtr partner = initial->progenitor()->partner(); Lorentz5Momentum ppartner[2]; if(partner) ppartner[0]=partner->momentum(); // get the n reference vector if(partner) { if(initial->progenitor()->showerKinematics()) { nvect = initial->progenitor()->showerBasis()->getBasis()[1]; } else { Lorentz5Momentum ppartner=initial->progenitor()->partner()->momentum(); if(gottaBoost) ppartner.boost(boosttorest,gammarest); nvect = Lorentz5Momentum( ZERO,0.5*initial->progenitor()->mass()* ppartner.vect().unit()); nvect.boost(-boosttorest,gammarest); } } // if ISR if(radiated[0]) { // reconstruct the decay jet reconstructDecayJet(initial->progenitor()); // momentum of decaying particle after ISR pjet=initial->progenitor()->momentum() -decay->incomingLines().begin()->second->momentum(); pjet.rescaleMass(); } // boost initial state jet and basis vector if needed if(gottaBoost) { pjet.boost(boosttorest,gammarest); nvect.boost(boosttorest,gammarest); ppartner[0].boost(boosttorest,gammarest); } // loop over the final-state particles and do the reconstruction JetKinVect possiblepartners; JetKinVect jetKinematics; bool atLeastOnce = radiated[0]; LorentzRotation restboost(boosttorest,gammarest); Energy inmass(ZERO); for(unsigned int ix=0;ixprogenitor()->isFinalState()) { inmass=ShowerHardJets[ix]->progenitor()->mass(); continue; } // do the reconstruction JetKinStruct tempJetKin; tempJetKin.parent = ShowerHardJets[ix]->progenitor(); if(ShowerHardJets.size()==2) { Lorentz5Momentum dum=ShowerHardJets[ix]->progenitor()->momentum(); dum.setMass(inmass); dum.rescaleRho(); tempJetKin.parent->set5Momentum(dum); } tempJetKin.p = ShowerHardJets[ix]->progenitor()->momentum(); if(gottaBoost) tempJetKin.p.boost(boosttorest,gammarest); _progenitor=tempJetKin.parent; if(ShowerHardJets[ix]->reconstructed()==ShowerProgenitor::notReconstructed) { atLeastOnce |= reconstructTimeLikeJet(tempJetKin.parent); ShowerHardJets[ix]->reconstructed(ShowerProgenitor::done); } if(gottaBoost) deepTransform(tempJetKin.parent,restboost); tempJetKin.q = ShowerHardJets[ix]->progenitor()->momentum(); jetKinematics.push_back(tempJetKin); } if(partner) ppartner[1]=partner->momentum(); // calculate the rescaling parameters double k1,k2; Lorentz5Momentum qt; if(!solveDecayKFactor(initial->progenitor()->mass(),nvect,pjet, jetKinematics,partner,ppartner,k1,k2,qt)) { for(map >::const_iterator bit=_boosts.begin();bit!=_boosts.end();++bit) { for(vector::const_reverse_iterator rit=bit->second.rbegin();rit!=bit->second.rend();++rit) { LorentzRotation rot = rit->inverse(); bit->first->transform(rot); } } _boosts.clear(); for(map >::const_iterator bit=_treeBoosts.begin();bit!=_treeBoosts.end();++bit) { for(vector::const_reverse_iterator rit=bit->second.rbegin();rit!=bit->second.rend();++rit) { LorentzRotation rot = rit->inverse(); bit->first->transform(rot,false); } } _treeBoosts.clear(); _currentTree = tShowerTreePtr(); return false; } // apply boosts and rescalings to final-state jets for(JetKinVect::iterator it = jetKinematics.begin(); it != jetKinematics.end(); ++it) { LorentzRotation Trafo = LorentzRotation(); if(it->parent!=partner) { // boost for rescaling if(atLeastOnce) { map >::const_iterator tit; for(tit = _currentTree->treelinks().begin(); tit != _currentTree->treelinks().end();++tit) { if(tit->second.first && tit->second.second==it->parent) break; } if(it->parent->children().empty()&&!it->parent->spinInfo() && tit==_currentTree->treelinks().end()) { Lorentz5Momentum pnew(k2*it->p.vect(), sqrt(sqr(k2*it->p.vect().mag())+it->q.mass2()), it->q.mass()); it->parent->set5Momentum(pnew); } else { // rescaling boost can't ever work in this case if(k2<0. && it->q.mass()==ZERO) throw KinematicsReconstructionVeto(); Trafo = solveBoost(k2, it->q, it->p); } } if(gottaBoost) Trafo.boost(-boosttorest,gammarest); if(atLeastOnce || gottaBoost) deepTransform(it->parent,Trafo); } else { Lorentz5Momentum pnew=ppartner[0]; pnew *=k1; pnew-=qt; pnew.setMass(ppartner[1].mass()); pnew.rescaleEnergy(); LorentzRotation Trafo=solveBoost(1.,ppartner[1],pnew); if(gottaBoost) Trafo.boost(-boosttorest,gammarest); deepTransform(partner,Trafo); } } } catch(KinematicsReconstructionVeto) { for(map >::const_iterator bit=_boosts.begin();bit!=_boosts.end();++bit) { for(vector::const_reverse_iterator rit=bit->second.rbegin();rit!=bit->second.rend();++rit) { LorentzRotation rot = rit->inverse(); bit->first->transform(rot); } } _boosts.clear(); for(map >::const_iterator bit=_treeBoosts.begin();bit!=_treeBoosts.end();++bit) { for(vector::const_reverse_iterator rit=bit->second.rbegin();rit!=bit->second.rend();++rit) { LorentzRotation rot = rit->inverse(); bit->first->transform(rot,false); } } _treeBoosts.clear(); _currentTree = tShowerTreePtr(); return false; } catch (Exception & ex) { _currentTree = tShowerTreePtr(); _boosts.clear(); _treeBoosts.clear(); throw ex; } _boosts.clear(); _treeBoosts.clear(); _currentTree = tShowerTreePtr(); return true; } bool KinematicsReconstructor:: reconstructDecayJet( const tShowerParticlePtr p) const { if(p->children().empty()) return false; tShowerParticlePtr child; // if branching reconstruct time-like child child = dynamic_ptr_cast(p->children()[1]); if(child) { _progenitor=child; reconstructTimeLikeJet(child); // calculate the momentum of the particle Lorentz5Momentum pnew=p->momentum()-child->momentum(); pnew.rescaleMass(); p->children()[0]->set5Momentum(pnew); child=dynamic_ptr_cast(p->children()[0]); reconstructDecayJet(child); return true; } return false; } bool KinematicsReconstructor:: solveDecayKFactor(Energy mb, const Lorentz5Momentum & n, const Lorentz5Momentum & pjet, const JetKinVect & jetKinematics, ShowerParticlePtr partner, Lorentz5Momentum ppartner[2], double & k1, double & k2, Lorentz5Momentum & qt) const { Energy2 pjn = partner ? pjet.vect()*n.vect() : ZERO; Energy2 pcn = partner ? ppartner[0].vect()*n.vect() : 1.*MeV2; Energy2 nmag = n.vect().mag2(); Lorentz5Momentum pn = partner ? (pjn/nmag)*n : Lorentz5Momentum(); qt=pjet-pn; qt.setE(ZERO); Energy2 pt2=qt.vect().mag2(); Energy Ejet = pjet.e(); // magnitudes of the momenta for fast access vector pmag; Energy total(Ejet); for(unsigned int ix=0;ixmb) return false; Energy2 pcmag=ppartner[0].vect().mag2(); // used newton-raphson to get the rescaling static const Energy eps=1e-8*GeV; long double d1(1.),d2(1.); Energy roots, ea, ec, ds; unsigned int ix=0; do { ++ix; d2 = d1 + pjn/pcn; roots = Ejet; ds = ZERO; for(unsigned int iy=0;iyeps && ix<100); k1=d1; k2=d2; // return true if N-R succeed, otherwise false return ix<100; } bool KinematicsReconstructor:: deconstructDecayJets(HardTreePtr decay,ShowerInteraction) const { // extract the momenta of the particles vector pin; vector pout; // on-shell masses of the decay products vector mon; Energy mbar(-GeV); // the hard branchings of the particles set::iterator cit; set branchings=decay->branchings(); // properties of the incoming particle bool ISR = false; HardBranchingPtr initial; Lorentz5Momentum qisr; // find the incoming particle, both before and after // any ISR for(cit=branchings.begin();cit!=branchings.end();++cit){ if((*cit)->status()==HardBranching::Incoming|| (*cit)->status()==HardBranching::Decay) { // search back up isr if needed HardBranchingPtr branch = *cit; while(branch->parent()) branch=branch->parent(); initial=branch; // momentum or original parent pin.push_back(branch->branchingParticle()->momentum()); // ISR? ISR = !branch->branchingParticle()->children().empty(); // ISR momentum qisr = pin.back()-(**cit).branchingParticle()->momentum(); qisr.rescaleMass(); } } assert(pin.size()==1); // compute boost to rest frame Boost boostv=-pin[0].boostVector(); // partner for ISR ShowerParticlePtr partner; Lorentz5Momentum ppartner; if(initial->branchingParticle()->partner()) { partner=initial->branchingParticle()->partner(); ppartner=partner->momentum(); } // momentum of the decay products for(cit=branchings.begin();cit!=branchings.end();++cit) { if((*cit)->status()!=HardBranching::Outgoing) continue; // find the mass of the particle // including special treatment for off-shell resonances // to preserve off-shell mass Energy mass; if(!(**cit).branchingParticle()->dataPtr()->stable()) { HardBranchingPtr branch=*cit; while(!branch->children().empty()) { for(unsigned int ix=0;ixchildren().size();++ix) { if(branch->children()[ix]->branchingParticle()->id()== (**cit).branchingParticle()->id()) { branch = branch->children()[ix]; continue; } } }; mass = branch->branchingParticle()->mass(); } else { mass = (**cit).branchingParticle()->dataPtr()->mass(); } // if not evolution partner of decaying particle if((*cit)->branchingParticle()!=partner) { pout.push_back((*cit)->branchingParticle()->momentum()); mon.push_back(mass); } // evolution partner of decaying particle else { mbar = mass; } } // boost all the momenta to the rest frame of the decaying particle for(unsigned int ix=0;ixbranchingParticle()->partner()) { ppartner.boost(boostv); qisr.boost(boostv); } // compute the rescaling factors double k1,k2; if(!ISR) { if(partner) { pout.push_back(ppartner); mon.push_back(mbar); } k1=k2=inverseRescalingFactor(pout,mon,pin[0].mass()); if(partner) { pout.pop_back(); mon.pop_back(); } } else { if(!inverseDecayRescalingFactor(pout,mon,pin[0].mass(), ppartner,mbar,k1,k2)) return false; } // now calculate the p reference vectors unsigned int ifinal=0; for(cit=branchings.begin();cit!=branchings.end();++cit) { if((**cit).status()!=HardBranching::Outgoing) continue; // for partners other than colour partner of decaying particle if((*cit)->branchingParticle()!=partner) { Lorentz5Momentum pvect = (*cit)->branchingParticle()->momentum(); pvect.boost(boostv); pvect /= k1; pvect.setMass(mon[ifinal]); ++ifinal; pvect.rescaleEnergy(); pvect.boost(-boostv); (*cit)->pVector(pvect); (*cit)->showerMomentum(pvect); } // for colour partner of decaying particle else { Lorentz5Momentum pvect = (*cit)->branchingParticle()->momentum(); pvect.boost(boostv); Lorentz5Momentum qtotal; for(unsigned int ix=0;ixpVector(pvect); (*cit)->showerMomentum(pvect); } } // For initial-state if needed if(initial) { tShowerParticlePtr newPartner=initial->branchingParticle()->partner(); if(newPartner) { tHardBranchingPtr branch; for( set::iterator clt = branchings.begin(); clt != branchings.end(); ++clt ) { if((**clt).branchingParticle()==newPartner) { initial->colourPartner(*clt); branch=*clt; break; } } Lorentz5Momentum pvect = initial->branchingParticle()->momentum(); initial->pVector(pvect); Lorentz5Momentum ptemp = branch->pVector(); ptemp.boost(boostv); Lorentz5Momentum nvect = Lorentz5Momentum( ZERO, 0.5*initial->branchingParticle()->mass()* ptemp.vect().unit()); nvect.boost(-boostv); initial->nVector(nvect); } } // calculate the reference vectors, then for outgoing particles for(cit=branchings.begin();cit!=branchings.end();++cit){ if((**cit).status()!=HardBranching::Outgoing) continue; // find the partner branchings tShowerParticlePtr newPartner=(*cit)->branchingParticle()->partner(); if(!newPartner) continue; tHardBranchingPtr branch; for( set::iterator clt = branchings.begin(); clt != branchings.end(); ++clt ) { if(cit==clt) continue; if((**clt).branchingParticle()==newPartner) { (**cit).colourPartner(*clt); branch=*clt; break; } } if((**decay->incoming().begin()).branchingParticle()==newPartner) { (**cit).colourPartner(*decay->incoming().begin()); branch = *decay->incoming().begin(); } // final-state colour partner if(branch->status()==HardBranching::Outgoing) { Boost boost=((*cit)->pVector()+branch->pVector()).findBoostToCM(); Lorentz5Momentum pcm = branch->pVector(); pcm.boost(boost); Lorentz5Momentum nvect = Lorentz5Momentum(ZERO,pcm.vect()); nvect.boost( -boost); (*cit)->nVector(nvect); } // initial-state colour partner else { Boost boost=branch->pVector().findBoostToCM(); Lorentz5Momentum pcm = (*cit)->pVector(); pcm.boost(boost); Lorentz5Momentum nvect = Lorentz5Momentum( ZERO, -pcm.vect()); nvect.boost( -boost); (*cit)->nVector(nvect); } } // now compute the new momenta // and calculate the shower variables for(cit=branchings.begin();cit!=branchings.end();++cit) { if((**cit).status()!=HardBranching::Outgoing) continue; LorentzRotation B=LorentzRotation(-boostv); LorentzRotation A=LorentzRotation(boostv),R; if((*cit)->branchingParticle()==partner) { Lorentz5Momentum qnew; Energy2 dot=(*cit)->pVector()*(*cit)->nVector(); double beta = 0.5*((*cit)->branchingParticle()->momentum().m2() -sqr((*cit)->pVector().mass()))/dot; qnew=(*cit)->pVector()+beta*(*cit)->nVector(); qnew.rescaleMass(); // compute the boost R=B*solveBoost(A*qnew,A*(*cit)->branchingParticle()->momentum())*A; } else { Lorentz5Momentum qnew; if((*cit)->branchingParticle()->partner()) { Energy2 dot=(*cit)->pVector()*(*cit)->nVector(); double beta = 0.5*((*cit)->branchingParticle()->momentum().m2() -sqr((*cit)->pVector().mass()))/dot; qnew=(*cit)->pVector()+beta*(*cit)->nVector(); qnew.rescaleMass(); } else { qnew = (*cit)->pVector(); } // compute the boost R=B*solveBoost(A*qnew,A*(*cit)->branchingParticle()->momentum())*A; } // reconstruct the momenta (*cit)->setMomenta(R,1.0,Lorentz5Momentum()); } if(initial) { initial->setMomenta(LorentzRotation(),1.0,Lorentz5Momentum()); } return true; } double KinematicsReconstructor:: inverseRescalingFactor(vector pout, vector mon, Energy roots) const { double lambda=1.; if(pout.size()==2) { double mu_q1(pout[0].m()/roots), mu_q2(pout[1].m()/roots); double mu_p1(mon[0]/roots) , mu_p2(mon[1]/roots); lambda = ((1.+mu_q1+mu_q2)*(1.-mu_q1-mu_q2)*(mu_q1-1.-mu_q2)*(mu_q2-1.-mu_q1))/ ((1.+mu_p1+mu_p2)*(1.-mu_p1-mu_p2)*(mu_p1-1.-mu_p2)*(mu_p2-1.-mu_p1)); if(lambda<0.) throw Exception() << "Rescaling factor is imaginary in KinematicsReconstructor::" << "inverseRescalingFactor lambda^2= " << lambda << Exception::eventerror; lambda = sqrt(lambda); } else { unsigned int ntry=0; // compute magnitudes once for speed vector pmag; for(unsigned int ix=0;ix root(pout.size()); do { // compute new energies Energy sum(ZERO); for(unsigned int ix=0;ix::const_iterator it=tree->branchings().begin(); it!=tree->branchings().end();++it) { if((**it).status()==HardBranching::Incoming) in .jets.push_back(*it); else out.jets.push_back(*it); } LorentzRotation toRest,fromRest; bool applyBoost(false); // do the initial-state reconstruction deconstructInitialInitialSystem(applyBoost,toRest,fromRest, tree,in.jets,type); // do the final-state reconstruction deconstructFinalStateSystem(toRest,fromRest,tree, out.jets,type); // only at this point that we can be sure all the reference vectors // are correct for(set::const_iterator it=tree->branchings().begin(); it!=tree->branchings().end();++it) { if((**it).status()==HardBranching::Incoming) continue; if((**it).branchingParticle()->coloured()) (**it).setMomenta(LorentzRotation(),1.,Lorentz5Momentum(),false); } for(set::const_iterator it=tree->incoming().begin(); it!=tree->incoming().end();++it) { (**it).setMomenta(LorentzRotation(),1.,Lorentz5Momentum(),false); } return true; } bool KinematicsReconstructor::deconstructHardJets(HardTreePtr tree, ShowerInteraction type) const { // inverse of old recon method if(_reconopt == 0) { return deconstructGeneralSystem(tree,type); } else if(_reconopt == 1) { return deconstructColourSinglets(tree,type); } else if(_reconopt == 2) { throw Exception() << "Inverse reconstruction is not currently supported for ReconstructionOption Colour2 " << "in KinematicsReconstructor::deconstructHardJets(). Please use one of the other options\n" << Exception::runerror; } else if(_reconopt == 3 || _reconopt == 4 ) { return deconstructColourPartner(tree,type); } else assert(false); } bool KinematicsReconstructor:: deconstructColourSinglets(HardTreePtr tree, ShowerInteraction type) const { // identify the colour singlet systems unsigned int nnun(0),nnii(0),nnif(0),nnf(0),nni(0); vector systems(identifySystems(tree->branchings(),nnun,nnii,nnif,nnf,nni)); // now decide what to do LorentzRotation toRest,fromRest; bool applyBoost(false); bool general(false); // initial-initial connection and final-state colour singlet systems // Drell-Yan type if(nnun==0&&nnii==1&&nnif==0&&nnf>0&&nni==0) { // reconstruct initial-initial system for(unsigned int ix=0;ix0&&nni==1)|| (nnif==2&& nni==0))) { for(unsigned int ix=0;ix0&&nni==2) { // only FS needed // but need to boost to rest frame if QED ISR Lorentz5Momentum ptotal; for(unsigned int ix=0;ixbranchingParticle()->momentum(); } toRest = LorentzRotation(ptotal.findBoostToCM()); fromRest = toRest; fromRest.invert(); if(type!=ShowerInteraction::QCD) { combineFinalState(systems); general=false; } } // general type else { general = true; } // final-state systems except for general recon if(!general) { for(unsigned int ix=0;ix::const_iterator it=tree->branchings().begin(); it!=tree->branchings().end();++it) { if((**it).status()==HardBranching::Incoming) continue; if((**it).branchingParticle()->coloured()) (**it).setMomenta(LorentzRotation(),1.,Lorentz5Momentum(),false); } for(set::const_iterator it=tree->incoming().begin(); it!=tree->incoming().end();++it) { (**it).setMomenta(LorentzRotation(),1.,Lorentz5Momentum(),false); } return true; } else { return deconstructGeneralSystem(tree,type); } return true; } bool KinematicsReconstructor:: deconstructColourPartner(HardTreePtr tree, ShowerInteraction type) const { Lorentz5Momentum ptotal; HardBranchingPtr emitter; ColourSingletShower incomingShower,outgoingShower; for(set::const_iterator it=tree->branchings().begin(); it!=tree->branchings().end();++it) { if((**it).status()==HardBranching::Incoming) { incomingShower.jets.push_back(*it); ptotal += (*it)->branchingParticle()->momentum(); // check for emitting particle if((**it).parent() ) { if(!emitter) emitter = *it; else throw Exception() << "Only one emitting particle allowed in " << "KinematicsReconstructor::deconstructColourPartner()" << Exception::runerror; } } else if ((**it).status()==HardBranching::Outgoing) { outgoingShower.jets.push_back(*it); // check for emitting particle if(!(**it).children().empty() ) { if(!emitter) emitter = *it; else throw Exception() << "Only one emitting particle allowed in " << "KinematicsReconstructor::deconstructColourPartner()" << Exception::runerror; } } } assert(emitter); assert(emitter->colourPartner()); ColourSingletShower system; system.jets.push_back(emitter); system.jets.push_back(emitter->colourPartner()); LorentzRotation toRest,fromRest; bool applyBoost(false); // identify the colour singlet system if(emitter->status() == HardBranching::Outgoing && emitter->colourPartner()->status() == HardBranching::Outgoing ) { system.type=F; // need to boost to rest frame if QED ISR if( !incomingShower.jets[0]->branchingParticle()->coloured() && !incomingShower.jets[1]->branchingParticle()->coloured() ) { Boost boost = ptotal.findBoostToCM(); toRest = LorentzRotation( boost); fromRest = LorentzRotation(-boost); } else findInitialBoost(ptotal,ptotal,toRest,fromRest); deconstructFinalStateSystem(toRest,fromRest,tree, system.jets,type); } else if (emitter->status() == HardBranching::Incoming && emitter->colourPartner()->status() == HardBranching::Incoming) { system.type=II; deconstructInitialInitialSystem(applyBoost,toRest,fromRest,tree,system.jets,type); // make sure the recoil gets applied deconstructFinalStateSystem(toRest,fromRest,tree, outgoingShower.jets,type); } else if ((emitter->status() == HardBranching::Outgoing && emitter->colourPartner()->status() == HardBranching::Incoming ) || (emitter->status() == HardBranching::Incoming && emitter->colourPartner()->status() == HardBranching::Outgoing)) { system.type=IF; // enusre incoming first if(system.jets[0]->status() == HardBranching::Outgoing) swap(system.jets[0],system.jets[1]); deconstructInitialFinalSystem(tree,system.jets,type); } else { throw Exception() << "Unknown type of system in " << "KinematicsReconstructor::deconstructColourPartner()" << Exception::runerror; } // only at this point that we can be sure all the reference vectors // are correct for(set::const_iterator it=tree->branchings().begin(); it!=tree->branchings().end();++it) { if((**it).status()==HardBranching::Incoming) continue; if((**it).branchingParticle()->coloured()) (**it).setMomenta(LorentzRotation(),1.,Lorentz5Momentum(),false); } for(set::const_iterator it=tree->incoming().begin(); it!=tree->incoming().end();++it) { (**it).setMomenta(LorentzRotation(),1.,Lorentz5Momentum(),false); } for(set::const_iterator it=tree->branchings().begin(); it!=tree->branchings().end();++it) { if((**it).status()!=HardBranching::Incoming) continue; if(*it==system.jets[0] || *it==system.jets[1]) continue; if((**it).branchingParticle()->momentum().z()>ZERO) { (**it).z((**it).branchingParticle()->momentum().plus()/(**it).beam()->momentum().plus()); } else { (**it).z((**it).branchingParticle()->momentum().minus()/(**it).beam()->momentum().minus()); } } return true; } void KinematicsReconstructor:: reconstructInitialFinalSystem(vector jets) const { Lorentz5Momentum pin[2],pout[2],pbeam; for(unsigned int ix=0;ixprogenitor()->isFinalState()) { pout[0] +=jets[ix]->progenitor()->momentum(); _progenitor = jets[ix]->progenitor(); if(jets[ix]->reconstructed()==ShowerProgenitor::notReconstructed) { reconstructTimeLikeJet(jets[ix]->progenitor()); jets[ix]->reconstructed(ShowerProgenitor::done); } } // initial-state parton else { pin[0] +=jets[ix]->progenitor()->momentum(); if(jets[ix]->progenitor()->showerKinematics()) { pbeam = jets[ix]->progenitor()->showerBasis()->getBasis()[0]; } else { if ( jets[ix]->original()->parents().empty() ) { pbeam = jets[ix]->progenitor()->momentum(); } else { pbeam = jets[ix]->original()->parents()[0]->momentum(); } } if(jets[ix]->reconstructed()==ShowerProgenitor::notReconstructed) { reconstructSpaceLikeJet(jets[ix]->progenitor()); jets[ix]->reconstructed(ShowerProgenitor::done); } assert(!jets[ix]->original()->parents().empty()); } } // add intrinsic pt if needed addIntrinsicPt(jets); // momenta after showering for(unsigned int ix=0;ixprogenitor()->isFinalState()) pout[1] += jets[ix]->progenitor()->momentum(); else pin[1] += jets[ix]->progenitor()->momentum(); } // work out the boost to the Breit frame Lorentz5Momentum pa = pout[0]-pin[0]; Axis axis(pa.vect().unit()); LorentzRotation rot; double sinth(sqrt(sqr(axis.x())+sqr(axis.y()))); if ( sinth > 1.e-9 ) rot.setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.)); rot.rotateX(Constants::pi); rot.boostZ( pa.e()/pa.vect().mag()); Lorentz5Momentum ptemp=rot*pbeam; Boost trans = -1./ptemp.e()*ptemp.vect(); trans.setZ(0.); if ( trans.mag2() - 1. >= 0. ) throw KinematicsReconstructionVeto(); rot.boost(trans); pa *=rot; // project and calculate rescaling // reference vectors Lorentz5Momentum n1(ZERO,ZERO,-pa.z(),-pa.z()); Lorentz5Momentum n2(ZERO,ZERO, pa.z(),-pa.z()); Energy2 n1n2 = n1*n2; // decompose the momenta Lorentz5Momentum qbp=rot*pin[1],qcp=rot*pout[1]; qbp.rescaleMass(); qcp.rescaleMass(); double a[2],b[2]; a[0] = n2*qbp/n1n2; b[0] = n1*qbp/n1n2; Lorentz5Momentum qperp = qbp-a[0]*n1-b[0]*n2; b[1] = 0.5; a[1] = 0.5*(qcp.m2()-qperp.m2())/n1n2/b[1]; double kb; if(a[0]!=0.) { double A(0.5*a[0]),B(b[0]*a[0]-a[1]*b[1]-0.25),C(-0.5*b[0]); if(sqr(B)-4.*A*C<0.) throw KinematicsReconstructionVeto(); kb = 0.5*(-B+sqrt(sqr(B)-4.*A*C))/A; } else { kb = 0.5*b[0]/(b[0]*a[0]-a[1]*b[1]-0.25); } // changed to improve stability if(kb==0.) throw KinematicsReconstructionVeto(); if ( a[1]>b[1] && abs(a[1]) < 1e-12 ) throw KinematicsReconstructionVeto(); if ( a[1]<=b[1] && abs(0.5+b[0]/kb) < 1e-12 ) throw KinematicsReconstructionVeto(); double kc = (a[1]>b[1]) ? (a[0]*kb-0.5)/a[1] : b[1]/(0.5+b[0]/kb); if(kc==0.) throw KinematicsReconstructionVeto(); Lorentz5Momentum pnew[2] = { a[0]*kb*n1+b[0]/kb*n2+qperp, a[1]*kc*n1+b[1]/kc*n2+qperp}; LorentzRotation rotinv=rot.inverse(); for(unsigned int ix=0;ixprogenitor()->isFinalState()) { deepTransform(jets[ix]->progenitor(),rot); deepTransform(jets[ix]->progenitor(),solveBoost(pnew[1],qcp)); Energy delta = jets[ix]->progenitor()->momentum().m()-jets[ix]->progenitor()->momentum().mass(); if ( abs(delta) > MeV ) throw KinematicsReconstructionVeto(); deepTransform(jets[ix]->progenitor(),rotinv); } else { tPPtr parent; boostChain(jets[ix]->progenitor(),rot,parent); boostChain(jets[ix]->progenitor(),solveBoostZ(pnew[0],qbp),parent); // check the first boost worked, and if not apply small correction to // fix energy/momentum conservation // this is a kludge but it reduces momentum non-conservation dramatically Lorentz5Momentum pdiff = pnew[0]-jets[ix]->progenitor()->momentum(); Energy2 delta = sqr(pdiff.x())+sqr(pdiff.y())+sqr(pdiff.z())+sqr(pdiff.t()); unsigned int ntry=0; while(delta>1e-6*GeV2 && ntry<5 ) { ntry +=1; boostChain(jets[ix]->progenitor(),solveBoostZ(pnew[0],jets[ix]->progenitor()->momentum()),parent); pdiff = pnew[0]-jets[ix]->progenitor()->momentum(); delta = sqr(pdiff.x())+sqr(pdiff.y())+sqr(pdiff.z())+sqr(pdiff.t()); } // apply test in breit-frame Lorentz5Momentum ptest1 = parent->momentum(); Lorentz5Momentum ptest2 = rot*pbeam; if(ptest1.z()/ptest2.z()<0. || ptest1.z()/ptest2.z()>1.) throw KinematicsReconstructionVeto(); boostChain(jets[ix]->progenitor(),rotinv,parent); } } } bool KinematicsReconstructor::addIntrinsicPt(vector jets) const { bool added=false; // add the intrinsic pt if needed for(unsigned int ix=0;ixprogenitor()->isFinalState()|| jets[ix]->hasEmitted()|| jets[ix]->reconstructed()==ShowerProgenitor::dontReconstruct) continue; if(_intrinsic.find(jets[ix])==_intrinsic.end()) continue; pair pt=_intrinsic[jets[ix]]; Energy etemp = jets[ix]->original()->parents()[0]->momentum().z(); Lorentz5Momentum p_basis(ZERO, ZERO, etemp, abs(etemp)), n_basis(ZERO, ZERO,-etemp, abs(etemp)); double alpha = jets[ix]->progenitor()->x(); double beta = 0.5*(sqr(jets[ix]->progenitor()->data().mass())+ sqr(pt.first))/alpha/(p_basis*n_basis); Lorentz5Momentum pnew=alpha*p_basis+beta*n_basis; pnew.setX(pt.first*cos(pt.second)); pnew.setY(pt.first*sin(pt.second)); pnew.rescaleMass(); jets[ix]->progenitor()->set5Momentum(pnew); added = true; } return added; } namespace { double defaultSolveBoostGamma(const double & betam,const Energy2 & kps, const Energy2 & qs, const Energy2 & Q2, const Energy & kp, const Energy & q, const Energy & qE) { if(betam<0.5) { return 1./sqrt(1.-sqr(betam)); } else { return ( kps+ qs + Q2)/ sqrt(2.*kps*qs + kps*Q2 + qs*Q2 + sqr(Q2) + 2.*q*qE*kp*sqrt(kps + Q2)); } } } LorentzRotation KinematicsReconstructor:: solveBoost(const double k, const Lorentz5Momentum & newq, const Lorentz5Momentum & oldp ) const { Energy q = newq.vect().mag(); Energy2 qs = sqr(q); Energy2 Q2 = newq.mass2(); Energy kp = k*(oldp.vect().mag()); Energy2 kps = sqr(kp); double betam = (q*newq.e() - kp*sqrt(kps + Q2))/(kps + qs + Q2); if ( abs(betam) - 1. >= 0. ) throw KinematicsReconstructionVeto(); Boost beta = -betam*(k/kp)*oldp.vect(); double gamma = 0.; if(Q2/sqr(oldp.e())>1e-4) { gamma = defaultSolveBoostGamma(betam,kps,qs,Q2,kp,q,newq.e()); } else { if(k>0) { gamma = 4.*kps*qs/sqr(kps +qs) + 2.*sqr(kps-qs)*Q2/pow<3,1>(kps +qs) - 0.25*( sqr(kps) + 14.*kps*qs + sqr(qs))*sqr(kps-qs)/(pow<4,1>(kps +qs)*kps*qs)*sqr(Q2); } else { gamma = 0.25*sqr(Q2)/(kps*qs)*(1. - 0.5*(kps+qs)/(kps*qs)*Q2); } if(gamma<=0.) throw KinematicsReconstructionVeto(); gamma = 1./sqrt(gamma); if(gamma>2.) gamma = defaultSolveBoostGamma(betam,kps,qs,Q2,kp,q,newq.e()); } // note that (k/kp)*oldp.vect() = oldp.vect()/oldp.vect().mag() but cheaper. ThreeVector ax = newq.vect().cross( oldp.vect() ); double delta; if (newq.x()*oldp.x()+newq.y()*oldp.y()+newq.z()*oldp.z()< 1e-16*GeV2) { throw KinematicsReconstructionVeto(); }else{ delta = newq.vect().angle( oldp.vect() ); } LorentzRotation R; using Constants::pi; Energy2 scale1 = sqr(newq.x())+ sqr(newq.y())+sqr(newq.z()); Energy2 scale2 = sqr(oldp.x())+ sqr(oldp.y())+sqr(oldp.z()); if ( ax.mag2()/scale1/scale2 > 1e-28 ) { R.rotate( delta, unitVector(ax) ).boost( beta , gamma ); } else if(abs(delta-pi)/pi < 0.001) { double phi=2.*pi*UseRandom::rnd(); Axis axis(cos(phi),sin(phi),0.); axis.rotateUz(newq.vect().unit()); R.rotate(delta,axis).boost( beta , gamma ); } else { R.boost( beta , gamma ); } return R; } LorentzRotation KinematicsReconstructor::solveBoost(const Lorentz5Momentum & q, const Lorentz5Momentum & p ) const { Energy modp = p.vect().mag(); Energy modq = q.vect().mag(); double betam = (p.e()*modp-q.e()*modq)/(sqr(modq)+sqr(modp)+p.mass2()); if ( abs(betam)-1. >= 0. ) throw KinematicsReconstructionVeto(); Boost beta = -betam*q.vect().unit(); ThreeVector ax = p.vect().cross( q.vect() ); double delta = p.vect().angle( q.vect() ); LorentzRotation R; using Constants::pi; if ( beta.mag2() - 1. >= 0. ) throw KinematicsReconstructionVeto(); if ( ax.mag2()/GeV2/MeV2 > 1e-16 ) { R.rotate( delta, unitVector(ax) ).boost( beta ); } else { R.boost( beta ); } return R; } LorentzRotation KinematicsReconstructor::solveBoostZ(const Lorentz5Momentum & q, const Lorentz5Momentum & p ) const { static const double eps = 1e-6; LorentzRotation R; double beta; Energy2 mt2 = p.mass()eps) { double erat = (q.t()+q.z())/(p.t()+p.z()); Energy2 den = mt2*(erat+1./erat); Energy2 num = (q.z()-p.z())*(q.t()+p.t()) + (p.z()+q.z())*(p.t()-q.t()); beta = num/den; if ( abs(beta) - 1. >= 0. ) throw KinematicsReconstructionVeto(); R.boostZ(beta); } else { double er = sqr(p.t()/q.t()); double x = ratio+0.125*(er+10.+1./er)*sqr(ratio); beta = -(p.t()-q.t())*(p.t()+q.t())/(sqr(p.t())+sqr(q.t()))*(1.+x); double gamma = (4.*sqr(p.t()*q.t()) +sqr(p.t()-q.t())*sqr(p.t()+q.t())* (-2.*x+sqr(x)))/sqr(sqr(p.t())+sqr(q.t())); if ( abs(beta) - 1. >= 0. ) throw KinematicsReconstructionVeto(); gamma = 1./sqrt(gamma); R.boost(0.,0.,beta,gamma); } Lorentz5Momentum ptest = R*p; if(ptest.z()/q.z() < 0. || ptest.t()/q.t() < 0. ) { throw KinematicsReconstructionVeto(); } return R; } void KinematicsReconstructor:: reconstructFinalStateSystem(bool applyBoost, const LorentzRotation & toRest, const LorentzRotation & fromRest, vector jets) const { LorentzRotation trans = applyBoost? toRest : LorentzRotation(); // special for case of individual particle if(jets.size()==1) { deepTransform(jets[0]->progenitor(),trans); deepTransform(jets[0]->progenitor(),fromRest); return; } bool radiated(false); // find the hard process centre-of-mass energy Lorentz5Momentum pcm; // check if radiated and calculate total momentum for(unsigned int ix=0;ixhasEmitted(); pcm += jets[ix]->progenitor()->momentum(); } if(applyBoost) pcm *= trans; // check if in CMF frame Boost beta_cm = pcm.findBoostToCM(); bool gottaBoost(false); if(beta_cm.mag() > 1e-12) { gottaBoost = true; trans.boost(beta_cm); } // collection of pointers to initial hard particle and jet momenta // for final boosts JetKinVect jetKinematics; vector::const_iterator cit; for(cit = jets.begin(); cit != jets.end(); cit++) { JetKinStruct tempJetKin; tempJetKin.parent = (*cit)->progenitor(); if(applyBoost || gottaBoost) { deepTransform(tempJetKin.parent,trans); } tempJetKin.p = (*cit)->progenitor()->momentum(); _progenitor=tempJetKin.parent; if((**cit).reconstructed()==ShowerProgenitor::notReconstructed) { radiated |= reconstructTimeLikeJet((*cit)->progenitor()); (**cit).reconstructed(ShowerProgenitor::done); } else { radiated |= !(*cit)->progenitor()->children().empty(); } tempJetKin.q = (*cit)->progenitor()->momentum(); jetKinematics.push_back(tempJetKin); } // default option rescale everything with the same factor if( _finalStateReconOption == 0 || jetKinematics.size() <= 2 ) { // find the rescaling factor double k = 0.0; if(radiated) { k = solveKfactor(pcm.m(), jetKinematics); // perform the rescaling and boosts for(JetKinVect::iterator it = jetKinematics.begin(); it != jetKinematics.end(); ++it) { LorentzRotation Trafo = solveBoost(k, it->q, it->p); deepTransform(it->parent,Trafo); } } } // different treatment of most off-shell else if ( _finalStateReconOption <= 4 ) { // sort the jets by virtuality std::sort(jetKinematics.begin(),jetKinematics.end(),JetOrdering()); // Bryan's procedures from FORTRAN if( _finalStateReconOption <=2 ) { // loop over the off-shell partons, _finalStateReconOption==1 only first ==2 all JetKinVect::const_iterator jend = _finalStateReconOption==1 ? jetKinematics.begin()+1 : jetKinematics.end(); for(JetKinVect::const_iterator jit=jetKinematics.begin(); jit!=jend;++jit) { // calculate the 4-momentum of the recoiling system Lorentz5Momentum psum; bool done = true; for(JetKinVect::const_iterator it=jetKinematics.begin();it!=jetKinematics.end();++it) { if(it==jit) { done = false; continue; } // first option put on-shell and sum 4-momenta if( _finalStateReconOption == 1 ) { LorentzRotation Trafo = solveBoost(1., it->q, it->p); deepTransform(it->parent,Trafo); psum += it->parent->momentum(); } // second option, sum momenta else { // already rescaled if(done) psum += it->parent->momentum(); // still needs to be rescaled else psum += it->p; } } // set the mass psum.rescaleMass(); // calculate the 3-momentum rescaling factor Energy2 s(pcm.m2()); Energy2 m1sq(jit->q.m2()),m2sq(psum.m2()); Energy4 num = sqr(s - m1sq - m2sq) - 4.*m1sq*m2sq; if(nump.vect().mag2()) ); // boost the off-shell parton LorentzRotation B1 = solveBoost(k, jit->q, jit->p); deepTransform(jit->parent,B1); // boost everything else to rescale LorentzRotation B2 = solveBoost(k, psum, psum); for(JetKinVect::iterator it=jetKinematics.begin();it!=jetKinematics.end();++it) { if(it==jit) continue; deepTransform(it->parent,B2); it->p *= B2; it->q *= B2; } } } // Peter's C++ procedures else { reconstructFinalFinalOffShell(jetKinematics,pcm.m2(), _finalStateReconOption == 4); } } else assert(false); // apply the final boosts if(gottaBoost || applyBoost) { LorentzRotation finalBoosts; if(gottaBoost) finalBoosts.boost(-beta_cm); if(applyBoost) finalBoosts.transform(fromRest); for(JetKinVect::iterator it = jetKinematics.begin(); it != jetKinematics.end(); ++it) { deepTransform(it->parent,finalBoosts); } } } void KinematicsReconstructor:: reconstructInitialInitialSystem(bool & applyBoost, LorentzRotation & toRest, LorentzRotation & fromRest, vector jets) const { bool radiated = false; Lorentz5Momentum pcm; // check whether particles radiated and calculate total momentum for( unsigned int ix = 0; ix < jets.size(); ++ix ) { radiated |= jets[ix]->hasEmitted(); pcm += jets[ix]->progenitor()->momentum(); if(jets[ix]->original()->parents().empty()) return; } pcm.rescaleMass(); // check if intrinsic pt to be added radiated |= !_intrinsic.empty(); // if no radiation return if(!radiated) { for(unsigned int ix=0;ixreconstructed()==ShowerProgenitor::notReconstructed) jets[ix]->reconstructed(ShowerProgenitor::done); } return; } // initial state shuffling applyBoost=false; vector p, pq, p_in; vector pts; for(unsigned int ix=0;ixprogenitor()->momentum()); // reconstruct the jet if(jets[ix]->reconstructed()==ShowerProgenitor::notReconstructed) { radiated |= reconstructSpaceLikeJet(jets[ix]->progenitor()); jets[ix]->reconstructed(ShowerProgenitor::done); } assert(!jets[ix]->original()->parents().empty()); Energy etemp = jets[ix]->original()->parents()[0]->momentum().z(); Lorentz5Momentum ptemp = Lorentz5Momentum(ZERO, ZERO, etemp, abs(etemp)); pq.push_back(ptemp); pts.push_back(jets[ix]->highestpT()); } // add the intrinsic pt if needed radiated |=addIntrinsicPt(jets); for(unsigned int ix=0;ixprogenitor()->momentum()); } double x1 = p_in[0].z()/pq[0].z(); double x2 = p_in[1].z()/pq[1].z(); vector beta=initialStateRescaling(x1,x2,p_in[0]+p_in[1],p,pq,pts); // if not need don't apply boosts if(!(radiated && p.size() == 2 && pq.size() == 2)) return; applyBoost=true; // apply the boosts Lorentz5Momentum newcmf; for(unsigned int ix=0;ixprogenitor(); Boost betaboost(0, 0, beta[ix]); tPPtr parent; boostChain(toBoost, LorentzRotation(0.,0.,beta[ix]),parent); if(parent->momentum().e()/pq[ix].e()>1.|| parent->momentum().z()/pq[ix].z()>1.) throw KinematicsReconstructionVeto(); newcmf+=toBoost->momentum(); } if(newcmf.m() jets, ShowerInteraction) const { assert(jets.size()==2); // put beam with +z first if(jets[0]->beam()->momentum().z() pin,pq; for(unsigned int ix=0;ixbranchingParticle()->momentum()); Energy etemp = jets[ix]->beam()->momentum().z(); pq.push_back(Lorentz5Momentum(ZERO, ZERO,etemp, abs(etemp))); } // calculate the rescaling double x[2]; Lorentz5Momentum pcm=pin[0]+pin[1]; assert(pcm.mass2()>ZERO); pcm.rescaleMass(); vector boost = inverseInitialStateRescaling(x[0],x[1],pcm,pin,pq); set::const_iterator cjt=tree->incoming().begin(); HardBranchingPtr incoming[2]; incoming[0] = *cjt; ++cjt; incoming[1] = *cjt; if((*tree->incoming().begin())->beam()->momentum().z()/pq[0].z()<0.) swap(incoming[0],incoming[1]); // apply the boost the the particles unsigned int iswap[2]={1,0}; for(unsigned int ix=0;ix<2;++ix) { LorentzRotation R(0.,0.,-boost[ix]); incoming[ix]->pVector(pq[ix]); incoming[ix]->nVector(pq[iswap[ix]]); incoming[ix]->setMomenta(R,1.,Lorentz5Momentum()); jets[ix]->showerMomentum(x[ix]*jets[ix]->pVector()); } // and calculate the boosts applyBoost=true; // do one boost if(_initialBoost==0) { toRest = LorentzRotation(-pcm.boostVector()); } else if(_initialBoost==1) { // first the transverse boost Energy pT = sqrt(sqr(pcm.x())+sqr(pcm.y())); double beta = -pT/pcm.t(); toRest=LorentzRotation(Boost(beta*pcm.x()/pT,beta*pcm.y()/pT,0.)); // the longitudinal beta = pcm.z()/sqrt(pcm.m2()+sqr(pcm.z())); toRest.boost(Boost(0.,0.,-beta)); } else assert(false); fromRest = LorentzRotation((jets[0]->showerMomentum()+ jets[1]->showerMomentum()).boostVector()); } void KinematicsReconstructor:: deconstructFinalStateSystem(const LorentzRotation & toRest, const LorentzRotation & fromRest, HardTreePtr tree, vector jets, ShowerInteraction type) const { LorentzRotation trans = toRest; if(jets.size()==1) { Lorentz5Momentum pnew = toRest*(jets[0]->branchingParticle()->momentum()); pnew *= fromRest; jets[0]-> original(pnew); jets[0]->showerMomentum(pnew); // find the colour partners ShowerParticleVector particles; vector ptemp; set::const_iterator cjt; for(cjt=tree->branchings().begin();cjt!=tree->branchings().end();++cjt) { ptemp.push_back((**cjt).branchingParticle()->momentum()); (**cjt).branchingParticle()->set5Momentum((**cjt).showerMomentum()); particles.push_back((**cjt).branchingParticle()); } dynamic_ptr_cast(ShowerHandler::currentHandler())->partnerFinder() ->setInitialEvolutionScales(particles,false,type,false); // calculate the reference vectors unsigned int iloc(0); set::iterator clt; for(cjt=tree->branchings().begin();cjt!=tree->branchings().end();++cjt) { // reset the momentum (**cjt).branchingParticle()->set5Momentum(ptemp[iloc]); ++iloc; // sort out the partners tShowerParticlePtr partner = (*cjt)->branchingParticle()->partner(); if(!partner) continue; for(clt=tree->branchings().begin();clt!=tree->branchings().end();++clt) { if((**clt).branchingParticle()==partner) { (**cjt).colourPartner(*clt); break; } } tHardBranchingPtr branch; for(clt=tree->branchings().begin();clt!=tree->branchings().end();++clt) { if(clt==cjt) continue; if((*clt)->branchingParticle()==partner) { branch=*clt; break; } } } return; } vector::iterator cit; vector pout; vector mon; Lorentz5Momentum pin; for(cit=jets.begin();cit!=jets.end();++cit) { pout.push_back((*cit)->branchingParticle()->momentum()); mon.push_back(findMass(*cit)); pin+=pout.back(); } // boost all the momenta to the rest frame of the decaying particle pin.rescaleMass(); pin *=trans; Boost beta_cm = pin.findBoostToCM(); bool gottaBoost(false); if(beta_cm.mag() > 1e-12) { gottaBoost = true; trans.boost(beta_cm); pin.boost(beta_cm); } for(unsigned int ix=0;ixbranchingParticle()->momentum(); pvect.transform(trans); pvect /= lambda; pvect.setMass(mon[ix]); pvect.rescaleEnergy(); if(gottaBoost) pvect.boost(-beta_cm); pvect.transform(fromRest); jets[ix]->pVector(pvect); jets[ix]->showerMomentum(pvect); } // find the colour partners ShowerParticleVector particles; vector ptemp; set::const_iterator cjt; for(cjt=tree->branchings().begin();cjt!=tree->branchings().end();++cjt) { ptemp.push_back((**cjt).branchingParticle()->momentum()); (**cjt).branchingParticle()->set5Momentum((**cjt).showerMomentum()); particles.push_back((**cjt).branchingParticle()); } dynamic_ptr_cast(ShowerHandler::currentHandler())->partnerFinder() ->setInitialEvolutionScales(particles,false,type,false); // calculate the reference vectors unsigned int iloc(0); set::iterator clt; for(cjt=tree->branchings().begin();cjt!=tree->branchings().end();++cjt) { // reset the momentum (**cjt).branchingParticle()->set5Momentum(ptemp[iloc]); ++iloc; } for(cjt=tree->branchings().begin();cjt!=tree->branchings().end();++cjt) { // sort out the partners tShowerParticlePtr partner = (*cjt)->branchingParticle()->partner(); if(!partner) continue; for(clt=tree->branchings().begin();clt!=tree->branchings().end();++clt) { if((**clt).branchingParticle()==partner) { (**cjt).colourPartner(*clt); break; } } tHardBranchingPtr branch; for(clt=tree->branchings().begin();clt!=tree->branchings().end();++clt) { if(clt==cjt) continue; if((*clt)->branchingParticle()==partner) { branch=*clt; break; } } // compute the reference vectors // both incoming, should all ready be done if((**cjt).status()==HardBranching::Incoming && (**clt).status()==HardBranching::Incoming) { continue; } // both outgoing else if((**cjt).status()!=HardBranching::Incoming&& branch->status()==HardBranching::Outgoing) { Boost boost=((*cjt)->pVector()+branch->pVector()).findBoostToCM(); Lorentz5Momentum pcm = branch->pVector(); pcm.boost(boost); Lorentz5Momentum nvect = Lorentz5Momentum(ZERO,pcm.vect()); nvect.boost( -boost); (**cjt).nVector(nvect); } else if((**cjt).status()==HardBranching::Incoming) { Lorentz5Momentum pa = -(**cjt).showerMomentum()+branch->showerMomentum(); Lorentz5Momentum pb = (**cjt).showerMomentum(); Axis axis(pa.vect().unit()); LorentzRotation rot; double sinth(sqrt(sqr(axis.x())+sqr(axis.y()))); rot.setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.)); rot.rotateX(Constants::pi); rot.boostZ( pa.e()/pa.vect().mag()); pb*=rot; Boost trans = -1./pb.e()*pb.vect(); trans.setZ(0.); rot.boost(trans); Energy scale=(**cjt).beam()->momentum().e(); Lorentz5Momentum pbasis(ZERO,(**cjt).beam()->momentum().vect().unit()*scale); Lorentz5Momentum pcm = rot*pbasis; rot.invert(); (**cjt).nVector(rot*Lorentz5Momentum(ZERO,-pcm.vect())); tHardBranchingPtr branch2 = *cjt;; while (branch2->parent()) { branch2=branch2->parent(); branch2->nVector(rot*Lorentz5Momentum(ZERO,-pcm.vect())); } } else if(branch->status()==HardBranching::Incoming) { (**cjt).nVector(Lorentz5Momentum(ZERO,branch->showerMomentum().vect())); } } // now compute the new momenta for(cjt=tree->branchings().begin();cjt!=tree->branchings().end();++cjt) { if(!(*cjt)->branchingParticle()->isFinalState()) continue; Lorentz5Momentum qnew; if((*cjt)->branchingParticle()->partner()) { Energy2 dot=(*cjt)->pVector()*(*cjt)->nVector(); double beta = 0.5*((*cjt)->branchingParticle()->momentum().m2() -sqr((*cjt)->pVector().mass()))/dot; qnew=(*cjt)->pVector()+beta*(*cjt)->nVector(); qnew.rescaleMass(); } else { qnew = (*cjt)->pVector(); } // qnew is the unshuffled momentum in the rest frame of the p basis vectors, // for the simple case Z->q qbar g this was checked against analytic formulae. // compute the boost LorentzRotation R=solveBoost(qnew, toRest*(*cjt)->branchingParticle()->momentum())*toRest; (*cjt)->setMomenta(R,1.0,Lorentz5Momentum()); } } Energy KinematicsReconstructor::momConsEq(double k, const Energy & root_s, const JetKinVect & jets) const { static const Energy2 eps=1e-8*GeV2; Energy dum = ZERO; for(JetKinVect::const_iterator it = jets.begin(); it != jets.end(); ++it) { Energy2 dum2 = (it->q).m2() + sqr(k)*(it->p).vect().mag2(); if(dum2 < ZERO) { if(dum2 < -eps) throw KinematicsReconstructionVeto(); dum2 = ZERO; } dum += sqrt(dum2); } return dum - root_s; } void KinematicsReconstructor::boostChain(tPPtr p, const LorentzRotation &bv, tPPtr & parent) const { if(!p->parents().empty()) boostChain(p->parents()[0], bv,parent); else parent=p; p->transform(bv); if(p->children().size()==2) { if(dynamic_ptr_cast(p->children()[1])) deepTransform(p->children()[1],bv); } } namespace { bool sortJets(ShowerProgenitorPtr j1, ShowerProgenitorPtr j2) { return j1->highestpT()>j2->highestpT(); } } void KinematicsReconstructor:: reconstructGeneralSystem(vector & ShowerHardJets) const { // find initial- and final-state systems ColourSingletSystem in,out; for(unsigned int ix=0;ixprogenitor()->isFinalState()) out.jets.push_back(ShowerHardJets[ix]); else in.jets.push_back(ShowerHardJets[ix]); } // reconstruct initial-initial system LorentzRotation toRest,fromRest; bool applyBoost(false); // reconstruct initial-initial system reconstructInitialInitialSystem(applyBoost,toRest,fromRest,in.jets); // reconstruct the final-state systems reconstructFinalStateSystem(applyBoost,toRest,fromRest,out.jets); } void KinematicsReconstructor:: reconstructFinalFirst(vector & ShowerHardJets) const { static const Energy2 minQ2 = 1e-4*GeV2; map used; for(unsigned int ix=0;ix outgoing; // first find any particles with final state partners for(unsigned int ix=0;ixprogenitor()->isFinalState()&& ShowerHardJets[ix]->progenitor()->partner()&& ShowerHardJets[ix]->progenitor()->partner()->isFinalState()) outgoing.insert(ShowerHardJets[ix]); } // then find the colour partners if(!outgoing.empty()) { set partners; for(set::const_iterator it=outgoing.begin();it!=outgoing.end();++it) { for(unsigned int ix=0;ixpartner()==ShowerHardJets[ix]->progenitor()) { partners.insert(ShowerHardJets[ix]); break; } } } outgoing.insert(partners.begin(),partners.end()); } // do the final-state reconstruction if needed if(!outgoing.empty()) { assert(outgoing.size()!=1); LorentzRotation toRest,fromRest; vector outgoingJets(outgoing.begin(),outgoing.end()); reconstructFinalStateSystem(false,toRest,fromRest,outgoingJets); } // Now do any initial-final systems which are needed vector IFSystems; // find the systems N.B. can have duplicates // find initial-state with FS partners or FS with IS partners for(unsigned int ix=0;ixprogenitor()->isFinalState()&& ShowerHardJets[ix]->progenitor()->partner()&& ShowerHardJets[ix]->progenitor()->partner()->isFinalState()) { IFSystems.push_back(ColourSingletSystem(IF,ShowerHardJets[ix])); } else if(ShowerHardJets[ix]->progenitor()->isFinalState()&& ShowerHardJets[ix]->progenitor()->partner()&& !ShowerHardJets[ix]->progenitor()->partner()->isFinalState()) { IFSystems.push_back(ColourSingletSystem(IF,ShowerHardJets[ix])); } } // then add the partners for(unsigned int is=0;isprogenitor()->partner()==ShowerHardJets[ix]->progenitor()) { IFSystems[is].jets.push_back(ShowerHardJets[ix]); } } // ensure incoming first if(IFSystems[is].jets[0]->progenitor()->isFinalState()) swap(IFSystems[is].jets[0],IFSystems[is].jets[1]); } if(!IFSystems.empty()) { unsigned int istart = UseRandom::irnd(IFSystems.size()); unsigned int istop=IFSystems.size(); for(unsigned int is=istart;is<=istop;++is) { if(is==IFSystems.size()) { if(istart!=0) { istop = istart-1; is=0; } else break; } // skip duplicates if(used[IFSystems[is].jets[0]] && used[IFSystems[is].jets[1]] ) continue; if(IFSystems[is].jets[0]->original()&&IFSystems[is].jets[0]->original()->parents().empty()) continue; Lorentz5Momentum psum; for(unsigned int ix=0;ixprogenitor()->isFinalState()) psum += IFSystems[is].jets[ix]->progenitor()->momentum(); else psum -= IFSystems[is].jets[ix]->progenitor()->momentum(); } if(-psum.m2()>minQ2) { reconstructInitialFinalSystem(IFSystems[is].jets); for(unsigned int ix=0;ixprogenitor()->isFinalState()) out.jets.push_back(ShowerHardJets[ix]); else in.jets.push_back(ShowerHardJets[ix]); } // reconstruct initial-initial system bool doRecon = false; for(unsigned int ix=0;ix & ShowerHardJets) const { static const Energy2 minQ2 = 1e-4*GeV2; // sort the vector by hardness of emission std::sort(ShowerHardJets.begin(),ShowerHardJets.end(),sortJets); // map between particles and progenitors for easy lookup map progenitorMap; for(unsigned int ix=0;ixprogenitor()] = ShowerHardJets[ix]; } // check that the IF systems can be reconstructed bool canReconstruct = true; for(unsigned int ix=0;ixprogenitor(); tShowerParticlePtr partner = progenitor->partner(); if(!partner) continue; else if((progenitor->isFinalState() && !partner->isFinalState()) || (!progenitor->isFinalState() && partner->isFinalState()) ) { vector jets(2); jets[0] = ShowerHardJets[ix]; jets[1] = progenitorMap[partner]; Lorentz5Momentum psum; for(unsigned int iy=0;iyprogenitor()->isFinalState()) psum += jets[iy]->progenitor()->momentum(); else psum -= jets[iy]->progenitor()->momentum(); } if(-psum.m2() used; for(unsigned int ix=0;ixreconstructed()==ShowerProgenitor::done) continue; // already reconstructed if(used[ShowerHardJets[ix]]) continue; // no partner continue tShowerParticlePtr progenitor = ShowerHardJets[ix]->progenitor(); tShowerParticlePtr partner = progenitor->partner(); if(!partner) { // check if there's a daughter tree which also needs boosting Lorentz5Momentum porig = progenitor->momentum(); map >::const_iterator tit; for(tit = _currentTree->treelinks().begin(); tit != _currentTree->treelinks().end();++tit) { // if there is, boost it if(tit->second.first && tit->second.second==progenitor) { Lorentz5Momentum pnew = tit->first->incomingLines().begin() ->first->progenitor()->momentum(); pnew *= tit->first->transform(); Lorentz5Momentum pdiff = porig-pnew; Energy2 test = sqr(pdiff.x()) + sqr(pdiff.y()) + sqr(pdiff.z()) + sqr(pdiff.t()); LorentzRotation rot; if(test>1e-6*GeV2) rot = solveBoost(porig,pnew); tit->first->transform(rot,false); _treeBoosts[tit->first].push_back(rot); } } ShowerHardJets[ix]->reconstructed(ShowerProgenitor::done); continue; } // do the reconstruction // final-final if(progenitor->isFinalState() && partner->isFinalState() ) { LorentzRotation toRest,fromRest; vector jets(2); jets[0] = ShowerHardJets[ix]; jets[1] = progenitorMap[partner]; if(_reconopt==4 && jets[1]->reconstructed()==ShowerProgenitor::notReconstructed) jets[1]->reconstructed(ShowerProgenitor::dontReconstruct); reconstructFinalStateSystem(false,toRest,fromRest,jets); if(_reconopt==4 && jets[1]->reconstructed()==ShowerProgenitor::dontReconstruct) jets[1]->reconstructed(ShowerProgenitor::notReconstructed); used[jets[0]] = true; if(_reconopt==3) used[jets[1]] = true; } // initial-final else if((progenitor->isFinalState() && !partner->isFinalState()) || (!progenitor->isFinalState() && partner->isFinalState()) ) { vector jets(2); jets[0] = ShowerHardJets[ix]; jets[1] = progenitorMap[partner]; if(jets[0]->progenitor()->isFinalState()) swap(jets[0],jets[1]); if(jets[0]->original()&&jets[0]->original()->parents().empty()) continue; Lorentz5Momentum psum; for(unsigned int iy=0;iyprogenitor()->isFinalState()) psum += jets[iy]->progenitor()->momentum(); else psum -= jets[iy]->progenitor()->momentum(); } if(_reconopt==4 && progenitorMap[partner]->reconstructed()==ShowerProgenitor::notReconstructed) progenitorMap[partner]->reconstructed(ShowerProgenitor::dontReconstruct); reconstructInitialFinalSystem(jets); if(_reconopt==4 && progenitorMap[partner]->reconstructed()==ShowerProgenitor::dontReconstruct) progenitorMap[partner]->reconstructed(ShowerProgenitor::notReconstructed); used[ShowerHardJets[ix]] = true; if(_reconopt==3) used[progenitorMap[partner]] = true; } // initial-initial else if(!progenitor->isFinalState() && !partner->isFinalState() ) { ColourSingletSystem in,out; in.jets.push_back(ShowerHardJets[ix]); in.jets.push_back(progenitorMap[partner]); for(unsigned int iy=0;iyprogenitor()->isFinalState()) out.jets.push_back(ShowerHardJets[iy]); } LorentzRotation toRest,fromRest; bool applyBoost(false); if(_reconopt==4 && in.jets[1]->reconstructed()==ShowerProgenitor::notReconstructed) in.jets[1]->reconstructed(ShowerProgenitor::dontReconstruct); reconstructInitialInitialSystem(applyBoost,toRest,fromRest,in.jets); if(_reconopt==4 && in.jets[1]->reconstructed()==ShowerProgenitor::dontReconstruct) in.jets[1]->reconstructed(ShowerProgenitor::notReconstructed); used[in.jets[0]] = true; if(_reconopt==3) used[in.jets[1]] = true; for(unsigned int iy=0;iyreconstructed()==ShowerProgenitor::notReconstructed) out.jets[iy]->reconstructed(ShowerProgenitor::dontReconstruct); } // reconstruct the final-state systems LorentzRotation finalBoosts; finalBoosts.transform( toRest); finalBoosts.transform(fromRest); for(unsigned int iy=0;iyprogenitor(),finalBoosts); } for(unsigned int iy=0;iyreconstructed()==ShowerProgenitor::dontReconstruct) out.jets[iy]->reconstructed(ShowerProgenitor::notReconstructed); } } } } bool KinematicsReconstructor:: inverseDecayRescalingFactor(vector pout, vector mon,Energy roots, Lorentz5Momentum ppartner, Energy mbar, double & k1, double & k2) const { ThreeVector qtotal; vector pmag; for(unsigned int ix=0;ix1e10) return false; } while (abs(numer)>eps&&itry<100); k1 = abs(k1); k2 = a*k1; return itry<100; } void KinematicsReconstructor:: deconstructInitialFinalSystem(HardTreePtr tree,vector jets, ShowerInteraction type) const { HardBranchingPtr incoming; Lorentz5Momentum pin[2],pout[2],pbeam; HardBranchingPtr initial; Energy mc(ZERO); for(unsigned int ix=0;ixstatus()==HardBranching::Outgoing) { pout[0] += jets[ix]->branchingParticle()->momentum(); mc = jets[ix]->branchingParticle()->thePEGBase() ? jets[ix]->branchingParticle()->thePEGBase()->mass() : jets[ix]->branchingParticle()->dataPtr()->mass(); } // initial-state parton else { pin[0] += jets[ix]->branchingParticle()->momentum(); initial = jets[ix]; pbeam = jets[ix]->beam()->momentum(); Energy scale=pbeam.t(); pbeam = Lorentz5Momentum(ZERO,pbeam.vect().unit()*scale); incoming = jets[ix]; while(incoming->parent()) incoming = incoming->parent(); } } if(jets.size()>2) { pout[0].rescaleMass(); mc = pout[0].mass(); } // work out the boost to the Breit frame Lorentz5Momentum pa = pout[0]-pin[0]; Axis axis(pa.vect().unit()); LorentzRotation rot; double sinth(sqrt(sqr(axis.x())+sqr(axis.y()))); if(axis.perp2()>0.) { rot.setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.)); rot.rotateX(Constants::pi); rot.boostZ( pa.e()/pa.vect().mag()); } // transverse part Lorentz5Momentum paxis=rot*pbeam; Boost trans = -1./paxis.e()*paxis.vect(); trans.setZ(0.); rot.boost(trans); pa *= rot; // reference vectors Lorentz5Momentum n1(ZERO,ZERO,-pa.z(),-pa.z()); Lorentz5Momentum n2(ZERO,ZERO, pa.z(),-pa.z()); Energy2 n1n2 = n1*n2; // decompose the momenta Lorentz5Momentum qbp=rot*pin[0],qcp= rot*pout[0]; double a[2],b[2]; a[0] = n2*qbp/n1n2; b[0] = n1*qbp/n1n2; a[1] = n2*qcp/n1n2; b[1] = n1*qcp/n1n2; Lorentz5Momentum qperp = qbp-a[0]*n1-b[0]*n2; // before reshuffling Energy Q = abs(pa.z()); double c = sqr(mc/Q); Lorentz5Momentum pb(ZERO,ZERO,0.5*Q*(1.+c),0.5*Q*(1.+c)); Lorentz5Momentum pc(ZERO,ZERO,0.5*Q*(c-1.),0.5*Q*(1.+c)); double anew[2],bnew[2]; anew[0] = pb*n2/n1n2; bnew[0] = 0.5*(qbp.m2()-qperp.m2())/n1n2/anew[0]; bnew[1] = pc*n1/n1n2; anew[1] = 0.5*qcp.m2()/bnew[1]/n1n2; Lorentz5Momentum qnewb = (anew[0]*n1+bnew[0]*n2+qperp); Lorentz5Momentum qnewc = (anew[1]*n1+bnew[1]*n2); // initial-state boost LorentzRotation rotinv=rot.inverse(); LorentzRotation transb=rotinv*solveBoostZ(qnewb,qbp)*rot; // final-state boost LorentzRotation transc=rotinv*solveBoost(qnewc,qcp)*rot; // this will need changing for more than one outgoing particle // set the pvectors for(unsigned int ix=0;ixstatus()==HardBranching::Incoming) { jets[ix]->pVector(pbeam); jets[ix]->showerMomentum(rotinv*pb); incoming->pVector(jets[ix]->pVector()); } else { jets[ix]->pVector(rotinv*pc); jets[ix]->showerMomentum(jets[ix]->pVector()); } } // find the colour partners ShowerParticleVector particles; vector ptemp; set::const_iterator cjt; for(cjt=tree->branchings().begin();cjt!=tree->branchings().end();++cjt) { ptemp.push_back((**cjt).branchingParticle()->momentum()); (**cjt).branchingParticle()->set5Momentum((**cjt).showerMomentum()); particles.push_back((**cjt).branchingParticle()); } dynamic_ptr_cast(ShowerHandler::currentHandler())->partnerFinder() ->setInitialEvolutionScales(particles,false,type,false); unsigned int iloc(0); for(cjt=tree->branchings().begin();cjt!=tree->branchings().end();++cjt) { // reset the momentum (**cjt).branchingParticle()->set5Momentum(ptemp[iloc]); ++iloc; } for(vector::const_iterator cjt=jets.begin(); cjt!=jets.end();++cjt) { // sort out the partners tShowerParticlePtr partner = (*cjt)->branchingParticle()->partner(); if(!partner) continue; tHardBranchingPtr branch; for(set::const_iterator clt=tree->branchings().begin();clt!=tree->branchings().end();++clt) { if((**clt).branchingParticle()==partner) { (**cjt).colourPartner(*clt); branch=*clt; break; } } // compute the reference vectors // both incoming, should all ready be done if((**cjt).status()==HardBranching::Incoming && branch->status()==HardBranching::Incoming) { Energy etemp = (*cjt)->beam()->momentum().z(); Lorentz5Momentum nvect(ZERO, ZERO,-etemp, abs(etemp)); tHardBranchingPtr branch2 = *cjt; (**cjt).nVector(nvect); while (branch2->parent()) { branch2=branch2->parent(); branch2->nVector(nvect); } } // both outgoing else if((**cjt).status()==HardBranching::Outgoing&& branch->status()==HardBranching::Outgoing) { Boost boost=((*cjt)->pVector()+branch->pVector()).findBoostToCM(); Lorentz5Momentum pcm = branch->pVector(); pcm.boost(boost); Lorentz5Momentum nvect = Lorentz5Momentum(ZERO,pcm.vect()); nvect.boost( -boost); (**cjt).nVector(nvect); } else if((**cjt).status()==HardBranching::Incoming) { Lorentz5Momentum pa = -(**cjt).showerMomentum()+branch->showerMomentum(); Lorentz5Momentum pb = (**cjt).showerMomentum(); Axis axis(pa.vect().unit()); LorentzRotation rot; double sinth(sqrt(sqr(axis.x())+sqr(axis.y()))); if(axis.perp2()>1e-20) { rot.setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.)); rot.rotateX(Constants::pi); } if(abs(1.-pa.e()/pa.vect().mag())>1e-6) rot.boostZ( pa.e()/pa.vect().mag()); pb*=rot; Boost trans = -1./pb.e()*pb.vect(); trans.setZ(0.); rot.boost(trans); Energy scale=(**cjt).beam()->momentum().t(); Lorentz5Momentum pbasis(ZERO,(**cjt).beam()->momentum().vect().unit()*scale); Lorentz5Momentum pcm = rot*pbasis; rot.invert(); Lorentz5Momentum nvect = rot*Lorentz5Momentum(ZERO,-pcm.vect()); (**cjt).nVector(nvect); tHardBranchingPtr branch2 = *cjt; while (branch2->parent()) { branch2=branch2->parent(); branch2->nVector(nvect); } } else if(branch->status()==HardBranching::Incoming) { Lorentz5Momentum nvect=Lorentz5Momentum(ZERO,branch->showerMomentum().vect()); (**cjt).nVector(nvect); } } // now compute the new momenta for(vector::const_iterator cjt=jets.begin(); cjt!=jets.end();++cjt) { if((**cjt).status()==HardBranching::Outgoing) { (**cjt).setMomenta(transc,1.,Lorentz5Momentum()); } } incoming->setMomenta(transb,1.,Lorentz5Momentum()); } void KinematicsReconstructor::deepTransform(PPtr particle, const LorentzRotation & r, bool match, PPtr original) const { if(_boosts.find(particle)!=_boosts.end()) { _boosts[particle].push_back(r); } Lorentz5Momentum porig = particle->momentum(); if(!original) original = particle; for ( int i = 0, N = particle->children().size(); i < N; ++i ) { deepTransform(particle->children()[i],r, particle->children()[i]->id()==original->id()&&match,original); } particle->transform(r); // transform the p and n vectors ShowerParticlePtr sparticle = dynamic_ptr_cast(particle); if(sparticle && sparticle->showerBasis()) { sparticle->showerBasis()->transform(r); } if ( particle->next() ) deepTransform(particle->next(),r,match,original); if(!match) return; if(!particle->children().empty()) return; // force the mass shell if(particle->dataPtr()->stable()) { Lorentz5Momentum ptemp = particle->momentum(); ptemp.rescaleEnergy(); particle->set5Momentum(ptemp); } // check if there's a daughter tree which also needs boosting map >::const_iterator tit; for(tit = _currentTree->treelinks().begin(); tit != _currentTree->treelinks().end();++tit) { // if there is, boost it if(tit->second.first && tit->second.second==original) { Lorentz5Momentum pnew = tit->first->incomingLines().begin() ->first->progenitor()->momentum(); pnew *= tit->first->transform(); Lorentz5Momentum pdiff = porig-pnew; Energy2 test = sqr(pdiff.x()) + sqr(pdiff.y()) + sqr(pdiff.z()) + sqr(pdiff.t()); LorentzRotation rot; if(test>1e-6*GeV2) rot = solveBoost(porig,pnew); tit->first->transform(r*rot,false); _treeBoosts[tit->first].push_back(r*rot); } } } void KinematicsReconstructor::reconstructFinalFinalOffShell(JetKinVect orderedJets, Energy2 s, bool recursive) const { JetKinVect::iterator jit; jit = orderedJets.begin(); ++jit; // 4-momentum of recoiling system Lorentz5Momentum psum; for( ; jit!=orderedJets.end(); ++jit) psum += jit->p; psum.rescaleMass(); // calculate the 3-momentum rescaling factor Energy2 m1sq(orderedJets.begin()->q.m2()),m2sq(psum.m2()); Energy4 num = sqr(s - m1sq - m2sq) - 4.*m1sq*m2sq; if(nump.vect().mag2()) ); // boost the most off-shell LorentzRotation B1 = solveBoost(k, orderedJets.begin()->q, orderedJets.begin()->p); deepTransform(orderedJets.begin()->parent,B1); // boost everything else // first to rescale LorentzRotation B2 = solveBoost(k, psum, psum); // and then to rest frame of new system Lorentz5Momentum pnew = B2*psum; pnew.rescaleMass(); B2.transform(pnew.findBoostToCM()); // apply transform (calling routine ensures at least 3 elements) jit = orderedJets.begin(); ++jit; for(;jit!=orderedJets.end();++jit) { deepTransform(jit->parent,B2); jit->p *= B2; jit->q *= B2; } JetKinVect newJets(orderedJets.begin()+1,orderedJets.end()); // final reconstruction if(newJets.size()==2 || !recursive ) { // rescaling factor double k = solveKfactor(psum.m(), newJets); // rescale jets in the new CMF for(JetKinVect::iterator it = newJets.begin(); it != newJets.end(); ++it) { LorentzRotation Trafo = solveBoost(k, it->q, it->p); deepTransform(it->parent,Trafo); } } // recursive else { std::sort(newJets.begin(),newJets.end(),JetOrdering()); reconstructFinalFinalOffShell(newJets,psum.m2(),recursive); } // finally boost back from new CMF LorentzRotation back(-pnew.findBoostToCM()); for(JetKinVect::iterator it = newJets.begin(); it != newJets.end(); ++it) { deepTransform(it->parent,back); } } Energy KinematicsReconstructor::findMass(HardBranchingPtr branch) const { // KH - 230909 - If the particle has no children then it will // not have showered and so it should be "on-shell" so we can // get it's mass from it's momentum. This means that the // inverseRescalingFactor doesn't give any nans or do things // it shouldn't if it gets e.g. two Z bosons generated with // off-shell masses. This is for sure not the best solution. // PR 1/1/10 modification to previous soln // PR 28/8/14 change to procedure and factorize into a function if(branch->children().empty()) { return branch->branchingParticle()->mass(); } else if(!branch->children().empty() && !branch->branchingParticle()->dataPtr()->stable() ) { for(unsigned int ix=0;ixchildren().size();++ix) { if(branch->branchingParticle()->id()== branch->children()[ix]->branchingParticle()->id()) return findMass(branch->children()[ix]); } } return branch->branchingParticle()->dataPtr()->mass(); } vector KinematicsReconstructor::inverseInitialStateRescaling(double & x1, double & x2, const Lorentz5Momentum & pold, const vector & p, const vector & pq) const { // hadronic CMS Energy2 s = (pq[0] +pq[1] ).m2(); // partonic CMS Energy MDY = pold.m(); // find alpha, beta and pt Energy2 p12=pq[0]*pq[1]; double a[2],b[2]; Lorentz5Momentum pt[2]; for(unsigned int ix=0;ix<2;++ix) { a[ix] = p[ix]*pq[1]/p12; b [ix] = p[ix]*pq[0]/p12; pt[ix] = p[ix]-a[ix]*pq[0]-b[ix]*pq[1]; } // compute kappa // we always want to preserve the mass of the system double k1(1.),k2(1.); if(_initialStateReconOption==0) { double rap=pold.rapidity(); x2 = MDY/sqrt(s*exp(2.*rap)); x1 = sqr(MDY)/s/x2; k1=a[0]/x1; k2=b[1]/x2; } // longitudinal momentum else if(_initialStateReconOption==1) { double A = 1.; double C = -sqr(MDY)/s; double B = 2.*pold.z()/sqrt(s); if(abs(B)>1e-10) { double discrim = 1.-4.*A*C/sqr(B); if(discrim < 0.) throw KinematicsReconstructionVeto(); x1 = B>0. ? 0.5*B/A*(1.+sqrt(discrim)) : 0.5*B/A*(1.-sqrt(discrim)); } else { x1 = -C/A; if( x1 <= 0.) throw KinematicsReconstructionVeto(); x1 = sqrt(x1); } x2 = sqr(MDY)/s/x1; k1=a[0]/x1; k2=b[1]/x2; } // preserve mass and don't scale the softer system // to reproduce the dipole kinematics else if(_initialStateReconOption==2) { // in this case kp = k1 or k2 depending on who's the harder guy k1 = a[0]*b[1]*s/sqr(MDY); if ( pt[0].perp2() < pt[1].perp2() ) swap(k1,k2); x1 = a[0]/k1; x2 = b[1]/k2; } else assert(false); // decompose the momenta double anew[2] = {a[0]/k1,a[1]*k2}; double bnew[2] = {b[0]*k1,b[1]/k2}; vector boost(2); for(unsigned int ix=0;ix<2;++ix) { boost[ix] = getBeta(a [ix]+b [ix], a[ix] -b [ix], anew[ix]+bnew[ix], anew[ix]-bnew[ix]); } return boost; } vector KinematicsReconstructor::initialStateRescaling(double x1, double x2, const Lorentz5Momentum & pold, const vector & p, const vector & pq, const vector& highestpts) const { Energy2 S = (pq[0]+pq[1]).m2(); // find alphas and betas in terms of desired basis Energy2 p12 = pq[0]*pq[1]; double a[2] = {p[0]*pq[1]/p12,p[1]*pq[1]/p12}; double b[2] = {p[0]*pq[0]/p12,p[1]*pq[0]/p12}; Lorentz5Momentum p1p = p[0] - a[0]*pq[0] - b[0]*pq[1]; Lorentz5Momentum p2p = p[1] - a[1]*pq[0] - b[1]*pq[1]; // compute kappa // we always want to preserve the mass of the system Energy MDY = pold.m(); Energy2 A = a[0]*b[1]*S; Energy2 B = Energy2(sqr(MDY)) - (a[0]*b[0]+a[1]*b[1])*S - (p1p+p2p).m2(); Energy2 C = a[1]*b[0]*S; double rad = 1.-4.*A*C/sqr(B); if(rad < 0.) throw KinematicsReconstructionVeto(); double kp = B/(2.*A)*(1.+sqrt(rad)); // now compute k1 // conserve rapidity double k1(0.); double k2(0.); if(_initialStateReconOption==0) { rad = kp*(b[0]+kp*b[1])/(kp*a[0]+a[1]); rad *= pq[0].z()1e-10) { double discrim = 1.-4.*a2*c2/sqr(b2); if(discrim < 0.) throw KinematicsReconstructionVeto(); k1 = b2>0. ? 0.5*b2/a2*(1.+sqrt(discrim)) : 0.5*b2/a2*(1.-sqrt(discrim)); } else { k1 = -c2/a2; if( k1 <= 0.) throw KinematicsReconstructionVeto(); k1 = sqrt(k1); } k2 = kp/k1; } // preserve mass and don't scale the softer system // to reproduce the dipole kinematics else if(_initialStateReconOption==2) { // in this case kp = k1 or k2 depending on who's the harder guy k1 = kp; k2 = 1.; if ( highestpts[0] < highestpts[1] ) swap(k1,k2); } else assert(false); // calculate the boosts vector beta(2); beta[0] = getBeta((a[0]+b[0]), (a[0]-b[0]), (k1*a[0]+b[0]/k1), (k1*a[0]-b[0]/k1)); beta[1] = getBeta((a[1]+b[1]), (a[1]-b[1]), (a[1]/k2+k2*b[1]), (a[1]/k2-k2*b[1])); if (pq[0].z() > ZERO) { beta[0] = -beta[0]; beta[1] = -beta[1]; } return beta; } void KinematicsReconstructor:: reconstructColourSinglets(vector & ShowerHardJets, ShowerInteraction type) const { // identify and catagorize the colour singlet systems unsigned int nnun(0),nnii(0),nnif(0),nnf(0),nni(0); vector systems(identifySystems(set(ShowerHardJets.begin(),ShowerHardJets.end()), nnun,nnii,nnif,nnf,nni)); // now decide what to do // initial-initial connection and final-state colour singlet systems LorentzRotation toRest,fromRest; bool applyBoost(false),general(false); // Drell-Yan type if(nnun==0&&nnii==1&&nnif==0&&nnf>0&&nni==0) { // reconstruct initial-initial system for(unsigned int ix=0;ix0&&nni==1)|| (nnif==2&& nni==0))) { // check these systems can be reconstructed for(unsigned int ix=0;ixprogenitor()->isFinalState()) q += systems[ix].jets[iy]->progenitor()->momentum(); else q -= systems[ix].jets[iy]->progenitor()->momentum(); } q.rescaleMass(); // check above cut if(abs(q.m())>=_minQ) continue; if(nnif==1&&nni==1) { throw KinematicsReconstructionVeto(); } else { general = true; break; } } if(!general) { for(unsigned int ix=0;ix0&&nni==2) { general = type!=ShowerInteraction::QCD; } // general type else { general = true; } // final-state systems except for general recon if(!general) { for(unsigned int ix=0;ix namespace Herwig { using namespace ThePEG; /**\ingroup Shower * Exception class * used to communicate failure of kinematics * reconstruction. */ struct KinematicsReconstructionVeto {}; /** \ingroup Shower * A simple struct to store the information we need on the * showering */ struct JetKinStruct { /** * Parent particle of the jet */ tShowerParticlePtr parent; /** * Momentum of the particle before reconstruction */ Lorentz5Momentum p; /** * Momentum of the particle after reconstruction */ Lorentz5Momentum q; }; /** * typedef for a vector of JetKinStruct */ typedef vector JetKinVect; /** * Enum to identify types of colour singlet systems */ enum SystemType { UNDEFINED=-1, II, IF, F ,I }; /** * Struct to store colour singlets */ template struct ColourSinglet { typedef vector > VecType; ColourSinglet() : type(UNDEFINED) {} ColourSinglet(SystemType intype,Value inpart) : type(intype),jets(1,inpart) {} /** * The type of system */ SystemType type; /** * The jets in the system */ vector jets; }; /** * Struct to store a colour singlet system of particles */ typedef ColourSinglet ColourSingletSystem; /** * Struct to store a colour singlet shower */ typedef ColourSinglet ColourSingletShower; /** \ingroup Shower * * This class is responsible for the kinematical reconstruction * after each showering step, and also for the necessary Lorentz boosts * in order to preserve energy-momentum conservation in the overall collision, * and also the invariant mass and the rapidity of the hard subprocess system. * In the case of multi-step showering, there will be not unnecessary * kinematical reconstructions. * * There is also the option of taking a set of momenta for the particles * and inverting the reconstruction to give the evolution variables for the * shower. * * Notice: * - although we often use the term "jet" in either methods or variables names, * or in comments, which could appear applicable only for QCD showering, * there is indeed no "dynamics" represented in this class: only kinematics * is involved, as the name of this class remainds. Therefore it can be used * for any kind of showers (QCD-,QED-,EWK-,... bremsstrahlung). * * @see ShowerParticle * @see ShowerKinematics * @see \ref KinematicsReconstructorInterfaces "The interfaces" * defined for KinematicsReconstructor. */ class KinematicsReconstructor: public Interfaced { public: /** * Default constructor */ KinematicsReconstructor() : _reconopt(0), _initialBoost(0), _finalStateReconOption(0), _initialStateReconOption(0), _minQ(MeV) {}; /** * Methods to reconstruct the kinematics of a scattering or decay process */ //@{ /** * Given in input a vector of the particles which initiated the showers * the method does the reconstruction of such jets, * including the appropriate boosts (kinematics reshufflings) * needed to conserve the total energy-momentum of the collision * and preserving the invariant mass and the rapidity of the * hard subprocess system. */ virtual bool reconstructHardJets(ShowerTreePtr hard, const map > & pt, ShowerInteraction type, bool switchRecon) const; /** * Given in input a vector of the particles which initiated the showers * the method does the reconstruction of such jets, * including the appropriate boosts (kinematics reshufflings) * needed to conserve the total energy-momentum of the collision * and preserving the invariant mass and the rapidity of the * hard subprocess system. */ virtual bool reconstructDecayJets(ShowerTreePtr decay, ShowerInteraction type) const; //@} /** * Methods to invert the reconstruction of the shower for * a scattering or decay process and calculate * the variables used to generate the * shower given the particles produced. * This is needed for the CKKW and POWHEG approaches */ //@{ /** * Given the particles, with a history which we wish to interpret * as a shower reconstruct the variables used to generate the * shower */ virtual bool deconstructDecayJets(HardTreePtr,ShowerInteraction) const; /** * Given the particles, with a history which we wish to interpret * as a shower reconstruct the variables used to generate the shower * for a hard process */ virtual bool deconstructHardJets(HardTreePtr,ShowerInteraction) const; //@} public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** * Methods to reconstruct the kinematics of individual jets */ //@{ /** * Given the particle (ShowerParticle object) that * originates a forward (time-like) jet, this method reconstructs the kinematics * of the jet. That is, by starting from the final grand-children (which * originates directly or indirectly from particleJetParent, * and which don't have children), and moving "backwards" (in a physical * time picture), towards the particleJetParent, the * ShowerKinematics objects associated with the various particles, * which have been created during the showering, are now completed. * In particular, at the end, we get the mass of the jet, which is the * main information we want. * This methods returns false if there was no radiation or rescaling required */ virtual bool reconstructTimeLikeJet(const tShowerParticlePtr particleJetParent) const; /** * Exactly similar to the previous one, but for a space-like jet. * Also in this case we start from the final grand-children (which * are childless) of the particle which originates the jet, but in * this case we proceed "forward" (in the physical time picture) * towards the particleJetParent. * This methods returns false if there was no radiation or rescaling required */ bool reconstructSpaceLikeJet(const tShowerParticlePtr particleJetParent) const; /** * Exactly similar to the previous one, but for a decay jet * This methods returns false if there was no radiation or rescaling required */ bool reconstructDecayJet(const tShowerParticlePtr particleJetParent) const; //@} /** * Methods to perform the reconstruction of various types of colour * singlet systems */ //@{ /** * Perform the reconstruction of a system with one incoming and at least one * outgoing particle */ void reconstructInitialFinalSystem(vector) const; /** * Perform the reconstruction of a system with only final-state * particles */ void reconstructFinalStateSystem(bool applyBoost, const LorentzRotation & toRest, const LorentzRotation & fromRest, vector) const; /** * Reconstruction of a general coloured system */ void reconstructGeneralSystem(vector & ShowerHardJets) const; /** * Reconstruction of a general coloured system doing * final-final, then initial-final and then initial-initial */ void reconstructFinalFirst(vector & ShowerHardJets) const; /** * Reconstruction of a general coloured system doing * colour parners */ void reconstructColourPartner(vector & ShowerHardJets) const; /** * Reconstruction based on colour singlet systems */ void reconstructColourSinglets(vector & ShowerHardJets, ShowerInteraction type) const; /** * Perform the reconstruction of a system with only final-state * particles */ void reconstructInitialInitialSystem(bool & applyBoost, LorentzRotation & toRest, LorentzRotation & fromRest, vector) const; //@} /** * Methods to perform the inverse reconstruction of various types of * colour singlet systems */ //@{ /** * Perform the inverse reconstruction of a system with only final-state * particles */ void deconstructFinalStateSystem(const LorentzRotation & toRest, const LorentzRotation & fromRest, HardTreePtr, vector, ShowerInteraction) const; /** * Perform the inverse reconstruction of a system with only initial-state * particles */ void deconstructInitialInitialSystem(bool & applyBoost, LorentzRotation & toRest, LorentzRotation & fromRest, HardTreePtr, vector, ShowerInteraction ) const; /** * Perform the inverse reconstruction of a system with only initial-state * particles */ void deconstructInitialFinalSystem(HardTreePtr, vector, ShowerInteraction ) const; bool deconstructGeneralSystem(HardTreePtr, ShowerInteraction) const; bool deconstructColourSinglets(HardTreePtr, ShowerInteraction) const; bool deconstructColourPartner(HardTreePtr, ShowerInteraction) const; //@} /** * Recursively treat the most off-shell paricle seperately * for final-final reconstruction */ void reconstructFinalFinalOffShell(JetKinVect orderedJets, Energy2 s, bool recursive) const; /** * Various methods for the Lorentz transforms needed to do the * rescalings */ //@{ /** * Compute the boost to get from the the old momentum to the new */ LorentzRotation solveBoost(const double k, const Lorentz5Momentum & newq, const Lorentz5Momentum & oldp) const; /** * Compute the boost to get from the the old momentum to the new */ LorentzRotation solveBoost(const Lorentz5Momentum & newq, const Lorentz5Momentum & oldq) const; /** * Compute the boost to get from the the old momentum to the new */ LorentzRotation solveBoostZ(const Lorentz5Momentum & newq, const Lorentz5Momentum & oldq) const; /** * Recursively boost the initial-state shower * @param p The particle * @param bv The boost * @param parent The parent of the chain */ void boostChain(tPPtr p, const LorentzRotation & bv, tPPtr & parent) const; /** * Given a 5-momentum and a scale factor, the method returns the * Lorentz boost that transforms the 3-vector vec{momentum} ---> * k*vec{momentum}. The method returns the null boost in the case no * solution exists. This will only work in the case where the * outgoing jet-momenta are parallel to the momenta of the particles * leaving the hard subprocess. */ Boost solveBoostBeta( const double k, const Lorentz5Momentum & newq, const Lorentz5Momentum & oldp); /** * Compute boost parameter along z axis to get (Ep, any perp, qp) * from (E, same perp, q). */ double getBeta(const double E, const double q, const double Ep, const double qp) const {return (q*E-qp*Ep)/(sqr(qp)+sqr(E));} //@} /** * Methods to calculate the various scaling factors */ //@{ /** * Given a vector of 5-momenta of jets, where the 3-momenta are the initial * ones before showering and the masses are reconstructed after the showering, * this method returns the overall scaling factor for the 3-momenta of the * vector of particles, vec{P}_i -> k * vec{P}_i, such to preserve energy- * momentum conservation, i.e. after the rescaling the center of mass 5-momentum * is equal to the one specified in input, cmMomentum. * The method returns 0 if such factor cannot be found. * @param root_s Centre-of-mass energy * @param jets The jets */ double solveKfactor( const Energy & root_s, const JetKinVect & jets ) const; /** * Calculate the rescaling factors for the jets in a particle decay where * there was initial-state radiation * @param mb The mass of the decaying particle * @param n The reference vector for the initial state radiation * @param pjet The momentum of the initial-state jet * @param jetKinematics The JetKinStruct objects for the jets * @param partner The colour partner * @param ppartner The momentum of the colour partner of the decaying particle * before and after radiation * @param k1 The rescaling parameter for the partner * @param k2 The rescaling parameter for the outgoing singlet * @param qt The transverse momentum vector */ bool solveDecayKFactor(Energy mb, const Lorentz5Momentum & n, const Lorentz5Momentum & pjet, const JetKinVect & jetKinematics, ShowerParticlePtr partner, Lorentz5Momentum ppartner[2], double & k1, double & k2, Lorentz5Momentum & qt) const; /** * Compute the momentum rescaling factor needed to invert the shower * @param pout The momenta of the outgoing particles * @param mon The on-shell masses * @param roots The mass of the decaying particle */ double inverseRescalingFactor(vector pout, vector mon,Energy roots) const; /** * Compute the momentum rescaling factor needed to invert the shower * @param pout The momenta of the outgoing particles * @param mon The on-shell masses * @param roots The mass of the decaying particle * @param ppartner The momentum of the colour partner * @param mbar The mass of the decaying particle * @param k1 The first scaling factor * @param k2 The second scaling factor */ bool inverseDecayRescalingFactor(vector pout, vector mon,Energy roots, Lorentz5Momentum ppartner, Energy mbar, double & k1, double & k2) const; /** * Check the rescaling conserves momentum * @param k The rescaling * @param root_s The centre-of-mass energy * @param jets The jets */ Energy momConsEq(double k, const Energy & root_s, const JetKinVect & jets) const; void findInitialBoost(const Lorentz5Momentum & pold, const Lorentz5Momentum & pnew, LorentzRotation & toRest, LorentzRotation & fromRest) const; //@} /** * Find the colour partners of a particle to identify the colour singlet * systems for the reconstruction. */ template void findPartners(Value branch,set & done, const set & branchings, vector & jets) const; /** * Add the intrinsic \f$p_T\f$ to the system if needed */ bool addIntrinsicPt(vector) const; /** * Apply a transform to the particle and any child, including child ShowerTree * objects * @param particle The particle * @param r The Lorentz transformation * @param match Whether or not to look at children etc * @param original The original particle */ void deepTransform(PPtr particle,const LorentzRotation & r, bool match=true,PPtr original=PPtr()) const; /** * Find the mass of a particle in the hard branching */ Energy findMass(HardBranchingPtr) const; /** * Calculate the initial-state rescaling factors */ vector initialStateRescaling(double x1, double x2, const Lorentz5Momentum & pold, const vector & p, const vector & pq, const vector& highespts) const; /** * Calculate the inverse of the initial-state rescaling factor */ vector inverseInitialStateRescaling(double & x1, double & x2, const Lorentz5Momentum & pold, const vector & p, const vector & pq) const; /** * Find the colour singlet systems */ template typename ColourSinglet::VecType identifySystems(set jets, unsigned int & nnun,unsigned int & nnii, unsigned int & nnif,unsigned int & nnf, unsigned int & nni) const; /** * Combine final-state colour systems */ template void combineFinalState(vector > & systems) const; protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const {return new_ptr(*this);} /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const {return new_ptr(*this);} //@} protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ KinematicsReconstructor & operator=(const KinematicsReconstructor &) = delete; private: /** * Option for handling the reconstruction */ unsigned int _reconopt; /** * Option for the boost for initial-initial reconstruction */ unsigned int _initialBoost; /** * Option for the reconstruction of final state systems */ unsigned int _finalStateReconOption; /** * Option for the initial state reconstruction */ unsigned int _initialStateReconOption; /** * Minimum invariant mass for initial-final dipoles to allow the * reconstruction */ Energy _minQ; /** * The progenitor of the jet currently being reconstructed */ mutable tShowerParticlePtr _progenitor; /** * Storage of the intrinsic \f$p_T\f$ */ mutable map > _intrinsic; /** * Current ShowerTree */ mutable tShowerTreePtr _currentTree; /** * Particles which shouldn't have their masses rescaled as * vector for the interface */ PDVector _noRescaleVector; /** * Particles which shouldn't have their masses rescaled as * set for quick access */ set _noRescale; /** * Storage of the boosts applied to enable resetting after failure */ mutable map > _boosts; /** * Storage of the boosts applied to enable resetting after failure */ mutable map > _treeBoosts; }; } -#include "KinematicsReconstructor.tcc" #endif /* HERWIG_KinematicsReconstructor_H */