Page MenuHomeHEPForge
Files F8310277

No OneTemporary
Actions

View Options

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 <complex>
#include <cmath>
#include <iostream>
#include <cstdlib>
using namespace std;
namespace MG5_sm_COLOREA
{
void ixxxxx(double p[4], double fmass, int nhel, int nsf, complex<double> fi[6])
{
complex<double> chi[2];
double sf[2], sfomega[2], omega[2], pp, pp3, sqp0p3, sqm[2];
int ip, im, nh;
fi[0] = complex<double> (-p[0] * nsf, -p[3] * nsf);
fi[1] = complex<double> (-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<double> (sqrt(pp3 * 0.5/pp), 0);
if (pp3 == 0.0)
{
chi[1] = complex<double> (-nh, 0);
}
else
{
chi[1] = complex<double> (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<double> (sqp0p3, 0.0);
if (sqp0p3 == 0.0)
{
chi[1] = complex<double> (-nhel * sqrt(2.0 * p[0]), 0.0);
}
else
{
chi[1] = complex<double> (nh * p[1], p[2])/sqp0p3;
}
if (nh == 1)
{
fi[2] = complex<double> (0.0, 0.0);
fi[3] = complex<double> (0.0, 0.0);
fi[4] = chi[0];
fi[5] = chi[1];
}
else
{
fi[2] = chi[1];
fi[3] = chi[0];
fi[4] = complex<double> (0.0, 0.0);
fi[5] = complex<double> (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<double>
tc[18])
{
complex<double> 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<double> (p[0] * nst, p[3] * nst);
ft[5][0] = complex<double> (p[1] * nst, p[2] * nst);
// construct eps+
if(nhel >= 0)
{
if(pp == 0)
{
ep[0] = complex<double> (0, 0);
ep[1] = complex<double> (-sqh, 0);
ep[2] = complex<double> (0, nst * sqh);
ep[3] = complex<double> (0, 0);
}
else
{
ep[0] = complex<double> (0, 0);
ep[3] = complex<double> (pt/pp * sqh, 0);
if(pt != 0)
{
pzpt = p[3]/(pp * pt) * sqh;
ep[1] = complex<double> (-p[1] * pzpt, -nst * p[2]/pt * sqh);
ep[2] = complex<double> (-p[2] * pzpt, nst * p[1]/pt * sqh);
}
else
{
ep[1] = complex<double> (-sqh, 0);
ep[2] = complex<double> (0, nst * Sgn(sqh, p[3]));
}
}
}
// construct eps-
if(nhel <= 0)
{
if(pp == 0)
{
em[0] = complex<double> (0, 0);
em[1] = complex<double> (sqh, 0);
em[2] = complex<double> (0, nst * sqh);
em[3] = complex<double> (0, 0);
}
else
{
em[0] = complex<double> (0, 0);
em[3] = complex<double> (-pt/pp * sqh, 0);
if(pt != 0)
{
pzpt = -p[3]/(pp * pt) * sqh;
em[1] = complex<double> (-p[1] * pzpt, -nst * p[2]/pt * sqh);
em[2] = complex<double> (-p[2] * pzpt, nst * p[1]/pt * sqh);
}
else
{
em[1] = complex<double> (sqh, 0);
em[2] = complex<double> (0, nst * Sgn(sqh, p[3]));
}
}
}
// construct eps0
if(std::labs(nhel) <= 1)
{
if(pp == 0)
{
e0[0] = complex<double> (0, 0);
e0[1] = complex<double> (0, 0);
e0[2] = complex<double> (0, 0);
e0[3] = complex<double> (1, 0);
}
else
{
emp = p[0]/(tmass * pp);
e0[0] = complex<double> (pp/tmass, 0);
e0[3] = complex<double> (p[3] * emp, 0);
if(pt != 0)
{
e0[1] = complex<double> (p[1] * emp, 0);
e0[2] = complex<double> (p[2] * emp, 0);
}
else
{
e0[1] = complex<double> (0, 0);
e0[2] = complex<double> (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<double> 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<double> (p[0] * nsv, p[3] * nsv);
vc[1] = complex<double> (p[1] * nsv, p[2] * nsv);
if (vmass != 0.0)
{
hel0 = 1.0 - std::abs(hel);
if(pp == 0.0)
{
vc[2] = complex<double> (0.0, 0.0);
vc[3] = complex<double> (-hel * sqh, 0.0);
vc[4] = complex<double> (0.0, nsvahl * sqh);
vc[5] = complex<double> (hel0, 0.0);
}
else
{
emp = p[0]/(vmass * pp);
vc[2] = complex<double> (hel0 * pp/vmass, 0.0);
vc[5] = complex<double> (hel0 * p[3] * emp + hel * pt/pp * sqh, 0.0);
if (pt != 0.0)
{
pzpt = p[3]/(pp * pt) * sqh * hel;
vc[3] = complex<double> (hel0 * p[1] * emp - p[1] * pzpt, -nsvahl *
p[2]/pt * sqh);
vc[4] = complex<double> (hel0 * p[2] * emp - p[2] * pzpt, nsvahl *
p[1]/pt * sqh);
}
else
{
vc[3] = complex<double> (-hel * sqh, 0.0);
vc[4] = complex<double> (0.0, nsvahl * Sgn(sqh, p[3]));
}
}
}
else
{
pp = p[0];
pt = sqrt((p[1] * p[1]) + (p[2] * p[2]));
vc[2] = complex<double> (0.0, 0.0);
vc[5] = complex<double> (hel * pt/pp * sqh, 0.0);
if (pt != 0.0)
{
pzpt = p[3]/(pp * pt) * sqh * hel;
vc[3] = complex<double> (-p[1] * pzpt, -nsv * p[2]/pt * sqh);
vc[4] = complex<double> (-p[2] * pzpt, nsv * p[1]/pt * sqh);
}
else
{
vc[3] = complex<double> (-hel * sqh, 0.0);
vc[4] = complex<double> (0.0, nsv * Sgn(sqh, p[3]));
}
}
return;
}
void sxxxxx(double p[4], int nss, complex<double> sc[3])
{
sc[2] = complex<double> (1.00, 0.00);
sc[0] = complex<double> (p[0] * nss, p[3] * nss);
sc[1] = complex<double> (p[1] * nss, p[2] * nss);
return;
}
void oxxxxx(double p[4], double fmass, int nhel, int nsf, complex<double> fo[6])
{
complex<double> chi[2];
double sf[2], sfomeg[2], omega[2], pp, pp3, sqp0p3, sqm[2];
int nh, ip, im;
fo[0] = complex<double> (p[0] * nsf, p[3] * nsf);
fo[1] = complex<double> (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<double> (sqrt(pp3 * 0.5/pp), 0.00);
if (pp3 == 0.00)
{
chi[1] = complex<double> (-nh, 0.00);
}
else
{
chi[1] = complex<double> (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<double> (sqp0p3, 0.00);
if(sqp0p3 == 0.000)
{
chi[1] = complex<double> (-nhel, 0.00) * sqrt(2.0 * p[0]);
}
else
{
chi[1] = complex<double> (nh * p[1], -p[2])/sqp0p3;
}
if(nh == 1)
{
fo[2] = chi[0];
fo[3] = chi[1];
fo[4] = complex<double> (0.00, 0.00);
fo[5] = complex<double> (0.00, 0.00);
}
else
{
fo[2] = complex<double> (0.00, 0.00);
fo[3] = complex<double> (0.00, 0.00);
fo[4] = chi[1];
fo[5] = chi[0];
}
}
return;
}
void FFV1P0_3(std::complex<double> F1[], std::complex<double> F2[],
std::complex<double> COUP, double M3, double W3, std::complex<double> V3[])
{
static std::complex<double> cI = std::complex<double> (0., 1.);
double P3[4];
std::complex<double> 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<double> F1[], std::complex<double> V3[],
std::complex<double> COUP, double M2, double W2, std::complex<double> F2[])
{
static std::complex<double> cI = std::complex<double> (0., 1.);
double P2[4];
std::complex<double> 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<double> F1[], std::complex<double> V3[],
std::complex<double> COUP1, std::complex<double> COUP2, double M2, double
W2, std::complex<double> F2[])
{
- static std::complex<double> cI = std::complex<double> (0., 1.);
- std::complex<double> Ftmp[6];
- double P2[4];
+ std::complex<double> Ftmp[6];
std::complex<double> 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<double> F1[], std::complex<double> V3[],
std::complex<double> COUP, double M2, double W2, std::complex<double> F2[])
{
static std::complex<double> cI = std::complex<double> (0., 1.);
double P2[4];
std::complex<double> 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<double> F1[], std::complex<double> F2[],
std::complex<double> V3[], std::complex<double> COUP, std::complex<double>
& vertex)
{
static std::complex<double> cI = std::complex<double> (0., 1.);
std::complex<double> 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<double> F1[], std::complex<double> F2[],
std::complex<double> V3[], std::complex<double> COUP1, std::complex<double>
COUP2, std::complex<double> & vertex)
{
- static std::complex<double> cI = std::complex<double> (0., 1.);
std::complex<double> tmp;
FFV2_0(F1, F2, V3, COUP1, vertex);
FFV5_0(F1, F2, V3, COUP2, tmp);
vertex = vertex + tmp;
}
void FFV5_1(std::complex<double> F2[], std::complex<double> V3[],
std::complex<double> COUP, double M1, double W1, std::complex<double> F1[])
{
static std::complex<double> cI = std::complex<double> (0., 1.);
double P1[4];
std::complex<double> 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<double> F1[], std::complex<double> F2[],
std::complex<double> V3[], std::complex<double> COUP, std::complex<double>
& vertex)
{
static std::complex<double> cI = std::complex<double> (0., 1.);
std::complex<double> 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<double> V2[], std::complex<double> V3[],
std::complex<double> V4[], std::complex<double> COUP, double M1, double W1,
std::complex<double> V1[])
{
static std::complex<double> cI = std::complex<double> (0., 1.);
std::complex<double> TMP12;
std::complex<double> TMP11;
double P1[4];
std::complex<double> 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<double> V2[], std::complex<double> V3[],
std::complex<double> V4[], std::complex<double> COUP, double M1, double W1,
std::complex<double> V1[])
{
static std::complex<double> cI = std::complex<double> (0., 1.);
std::complex<double> TMP12;
double P1[4];
std::complex<double> TMP6;
std::complex<double> 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<double> V1[], std::complex<double> V2[],
std::complex<double> V3[], std::complex<double> COUP, std::complex<double>
& vertex)
{
static std::complex<double> cI = std::complex<double> (0., 1.);
std::complex<double> TMP2;
std::complex<double> TMP1;
double P1[4];
std::complex<double> TMP0;
double P2[4];
std::complex<double> TMP7;
double P3[4];
std::complex<double> TMP6;
std::complex<double> TMP5;
std::complex<double> TMP4;
std::complex<double> TMP3;
std::complex<double> 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<double> F1[], std::complex<double> F2[],
std::complex<double> COUP, double M3, double W3, std::complex<double> V3[])
{
static std::complex<double> cI = std::complex<double> (0., 1.);
std::complex<double> denom;
std::complex<double> 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<double> F1[], std::complex<double> F2[],
std::complex<double> COUP1, std::complex<double> COUP2, double M3, double
W3, std::complex<double> V3[])
{
- static std::complex<double> cI = std::complex<double> (0., 1.);
- std::complex<double> denom;
- double P3[4];
- double OM3;
+ std::complex<double> denom;
int i;
std::complex<double> 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<double> F1[], std::complex<double> V3[],
std::complex<double> COUP, double M2, double W2, std::complex<double> F2[])
{
static std::complex<double> cI = std::complex<double> (0., 1.);
double P2[4];
std::complex<double> 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<double> F2[], std::complex<double> V3[],
std::complex<double> COUP, double M1, double W1, std::complex<double> F1[])
{
static std::complex<double> cI = std::complex<double> (0., 1.);
double P1[4];
std::complex<double> 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<double> F2[], std::complex<double> V3[],
std::complex<double> COUP1, std::complex<double> COUP2, double M1, double
W1, std::complex<double> F1[])
{
- static std::complex<double> cI = std::complex<double> (0., 1.);
- double P1[4];
std::complex<double> denom;
int i;
std::complex<double> 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<double> F1[], std::complex<double> F2[],
std::complex<double> V3[], std::complex<double> COUP, std::complex<double>
& vertex)
{
static std::complex<double> cI = std::complex<double> (0., 1.);
std::complex<double> TMP15;
std::complex<double> 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<double> F2[], std::complex<double> V3[],
std::complex<double> COUP, double M1, double W1, std::complex<double> F1[])
{
static std::complex<double> cI = std::complex<double> (0., 1.);
double P1[4];
std::complex<double> 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<double> F1[], std::complex<double> F2[],
std::complex<double> COUP, double M3, double W3, std::complex<double> V3[])
{
static std::complex<double> cI = std::complex<double> (0., 1.);
std::complex<double> denom;
std::complex<double> TMP10;
double P3[4];
double OM3;
std::complex<double> 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<double> V2[], std::complex<double> V3[],
std::complex<double> V4[], std::complex<double> COUP, double M1, double W1,
std::complex<double> V1[])
{
static std::complex<double> cI = std::complex<double> (0., 1.);
std::complex<double> TMP11;
double P1[4];
std::complex<double> TMP6;
std::complex<double> 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<double> V2[], std::complex<double> V3[],
std::complex<double> COUP, double M1, double W1, std::complex<double> V1[])
{
static std::complex<double> cI = std::complex<double> (0., 1.);
std::complex<double> TMP2;
double P1[4];
std::complex<double> TMP0;
double P2[4];
double P3[4];
std::complex<double> TMP6;
std::complex<double> TMP5;
std::complex<double> TMP4;
std::complex<double> 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 <algorithm>
#include <iostream>
#include <fstream>
#include "read_slha_COLOREA.h"
void SLHABlock_COLOREA::set_entry(std::vector<int> 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<int> 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<int> 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<int> 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<int> 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<int> 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<int> 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<int> 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 <cassert>
#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;ix<beam->children().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<PPtr> original,copy;
for(unsigned int ix=0;ix<process->incoming().size();++ix) {
original.push_back(process->incoming()[ix].first);
copy .push_back(new_ptr(Particle(*original.back())));
}
for(unsigned int ix=0;ix<process->outgoing().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();ix<original.size();++ix) {
ShowerParticlePtr stemp= new_ptr(ShowerParticle(*copy[ix],itype,true));
fixColour(stemp);
_outgoingLines.insert(make_pair(new_ptr(ShowerProgenitor(original[ix],copy[ix],
stemp)),
stemp));
_forward.insert(stemp);
}
}
void ShowerTree::updateFinalStateShowerProduct(ShowerProgenitorPtr progenitor,
ShowerParticlePtr parent,
const ShowerParticleVector & children) {
assert(children.size()==2);
bool matches[2];
for(unsigned int ix=0;ix<2;++ix) {
matches[ix] = children[ix]->id()==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<ShowerProgenitorPtr,ShowerParticlePtr>::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<ShowerProgenitorPtr,tShowerParticlePtr>::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;ix<intermediates.size();++ix) {
init->abandonChild(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;ix<intermediates.size();++ix) {
init->abandonChild(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;ix<parents.size();++ix)
parents[ix]->abandonChild(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<ShowerParticlePtr>(*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<PPtr> original,copy;
for(unsigned int ix=0;ix<newProcess->incoming().size();++ix) {
original.push_back(newProcess->incoming()[ix].first);
copy .push_back(new_ptr(Particle(*original.back())));
}
for(unsigned int ix=0;ix<newProcess->outgoing().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<original.size();++ix) {
ShowerParticlePtr stemp= new_ptr(ShowerParticle(*copy[ix],2,true));
fixColour(stemp);
_outgoingLines.insert(make_pair(new_ptr(ShowerProgenitor(original[ix],copy[ix],
stemp)),
stemp));
_forward.insert(stemp);
}
}
void ShowerTree::insertDecay(StepPtr pstep,bool ISR, bool) {
assert(_incomingLines.size()==1);
colourLines().clear();
// find final particle from previous tree
PPtr final;
if(_parent&&!_parent->_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;ix<dpar.size();++ix) dpar[ix]->abandonChild(copy);
final->addChild(copy);
pstep->addDecayProduct(copy);
// just a copy does not move
copy->setLifeLength(Lorentz5Distance());
copy->setVertex(final->decayVertex());
// final-state radiation
map<ShowerProgenitorPtr,tShowerParticlePtr>::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;ix<parents.size();++ix)
{parents[ix]->abandonChild(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;ix<parents.size();++ix)
{parents[ix]->abandonChild(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<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cit;
map<ShowerProgenitorPtr, ShowerParticlePtr>::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;ix<children.size();++ix) orig->abandonChild(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;ix<children.size();++ix) orig->abandonChild(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<ShowerProgenitorPtr, ShowerParticlePtr>::const_iterator cit;
map<ShowerProgenitorPtr,tShowerParticlePtr>::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<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator mit;
map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::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<tShowerParticlePtr> 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<tShowerParticlePtr>::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<children.size(); ++ix) {
_forward.insert(children[ix]);
}
}
-void ShowerTree::addInitialStateBranching(ShowerParticlePtr oldParent,
- ShowerParticlePtr newParent,
+void ShowerTree::addInitialStateBranching(ShowerParticlePtr , ShowerParticlePtr ,
ShowerParticlePtr otherChild) {
_forward.insert(otherChild);
}
void ShowerTree::setParents() {
// set the parent tree of the children
map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator tit;
for(tit=_treelinks.begin();tit!=_treelinks.end();++tit)
tit->first->_parent=this;
}
vector<ShowerProgenitorPtr> ShowerTree::extractProgenitors() {
// extract the particles from the ShowerTree
map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator mit;
vector<ShowerProgenitorPtr> ShowerHardJets;
for(mit=incomingLines().begin();mit!=incomingLines().end();++mit)
ShowerHardJets.push_back((*mit).first);
map<ShowerProgenitorPtr,tShowerParticlePtr>::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<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator cit;
// incoming
for(cit=_incomingLines.begin();cit!=_incomingLines.end();++cit) {
cit->first->progenitor()->deepTransform(boost);
cit->first->copy()->deepTransform(boost);
}
// outgoing
map<ShowerProgenitorPtr,tShowerParticlePtr>::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<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator
tit=_treelinks.begin();tit!=_treelinks.end();++tit)
tit->first->transform(boost,applyNow);
}
void ShowerTree::applyTransforms() {
// now boost all the particles
map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator cit;
// incoming
for(cit=_incomingLines.begin();cit!=_incomingLines.end();++cit) {
cit->first->progenitor()->deepTransform(_transforms);
cit->first->copy()->deepTransform(_transforms);
}
// outgoing
map<ShowerProgenitorPtr,tShowerParticlePtr>::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<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator
tit=_treelinks.begin();tit!=_treelinks.end();++tit)
tit->first->applyTransforms();
_transforms = LorentzRotation();
}
void ShowerTree::clearTransforms() {
_transforms = LorentzRotation();
// // child trees
// for(map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::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<MultiColour>::pointer colour =
dynamic_ptr_cast<Ptr<MultiColour>::pointer>(part->colourInfo());
vector<tcColinePtr> lines = colour->colourLines();
for(unsigned int ix=0;ix<lines.size();++ix) {
ColinePtr line = const_ptr_cast<ColinePtr>(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<MultiColour>::pointer colour =
dynamic_ptr_cast<Ptr<MultiColour>::pointer>(part->colourInfo());
vector<tcColinePtr> lines = colour->antiColourLines();
for(unsigned int ix=0;ix<lines.size();++ix) {
ColinePtr line = const_ptr_cast<ColinePtr>(lines[ix]);
if(line) {
line->removeAntiColoured(part);
line->addAntiColoured(part);
}
}
}
}
}
vector<ShowerParticlePtr> ShowerTree::extractProgenitorParticles() {
vector<ShowerParticlePtr> particles;
// incoming particles
for(map<ShowerProgenitorPtr, ShowerParticlePtr>::const_iterator
cit=incomingLines().begin(); cit!=incomingLines().end();++cit)
particles.push_back(cit->first->progenitor());
// outgoing particles
for(map<ShowerProgenitorPtr,tShowerParticlePtr>::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<PerturbativeProcessPtr,ShowerTreePtr> 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<PerturbativeProcessPtr,ShowerTreePtr>::const_iterator it=treeMap.begin();
it!=treeMap.end();++it) {
// links to daughter trees
map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator mit;
for(mit=it->second->_outgoingLines.begin();
mit!=it->second->_outgoingLines.end();++mit) {
unsigned int iloc=0;
for(;iloc<it->first->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;ix<particle->children().size();++ix) {
output += sumMomenta(particle->children()[ix]);
}
return output;
}
}
void ShowerTree::checkMomenta() {
vector<Lorentz5Momentum> pin;
for(map<ShowerProgenitorPtr,ShowerParticlePtr>::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<Lorentz5Momentum> pout;
for(map<ShowerProgenitorPtr,tShowerParticlePtr>::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<double,double> x;
for(map<ShowerProgenitorPtr,ShowerParticlePtr>::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<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator it= _outgoingLines.begin();
it!=_outgoingLines.end();++it) {
output->bornOutgoing().push_back(it->first->progenitor());
}
return output;
}
void ShowerTree::setVetoes(const map<ShowerInteraction,Energy> & pTs,
unsigned int type) {
if(type==1||type==3) {
for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator it=_incomingLines.begin();
it!=_incomingLines.end();++it) {
for(map<ShowerInteraction,Energy>::const_iterator jt=pTs.begin();jt!=pTs.end();++jt)
it->first->maximumpT(jt->second,jt->first);
}
}
if(type==2||type==3) {
for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator it= _outgoingLines.begin();
it!=_outgoingLines.end();++it) {
for(map<ShowerInteraction,Energy>::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 <cassert>
+#include "KinematicsReconstructor.tcc"
using namespace Herwig;
DescribeClass<KinematicsReconstructor,Interfaced>
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<KinematicsReconstructor> 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<KinematicsReconstructor,unsigned int> 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<KinematicsReconstructor,Energy> 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<KinematicsReconstructor,ParticleData> 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<KinematicsReconstructor,unsigned int> 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<KinematicsReconstructor,unsigned int> 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<KinematicsReconstructor,unsigned int> 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<cPDPtr>(_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<ShowerParticlePtr>(*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<ShowerParticlePtr>
(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<tShowerProgenitorPtr,
pair<Energy,double> > & intrinsic,
ShowerInteraction type,
bool switchRecon) const {
_currentTree = hard;
_intrinsic=intrinsic;
// extract the particles from the ShowerTree
vector<ShowerProgenitorPtr> ShowerHardJets=hard->extractProgenitors();
for(unsigned int ix=0;ix<ShowerHardJets.size();++ix) {
_boosts[ShowerHardJets[ix]->progenitor()] = vector<LorentzRotation>();
}
for(map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator
tit = _currentTree->treelinks().begin();
tit != _currentTree->treelinks().end();++tit) {
_treeBoosts[tit->first] = vector<LorentzRotation>();
}
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<tPPtr,vector<LorentzRotation> >::const_iterator bit=_boosts.begin();bit!=_boosts.end();++bit) {
for(vector<LorentzRotation>::const_reverse_iterator rit=bit->second.rbegin();rit!=bit->second.rend();++rit) {
LorentzRotation rot = rit->inverse();
bit->first->transform(rot);
}
}
_boosts.clear();
for(map<tShowerTreePtr,vector<LorentzRotation> >::const_iterator bit=_treeBoosts.begin();bit!=_treeBoosts.end();++bit) {
for(vector<LorentzRotation>::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<ShowerProgenitorPtr,ShowerParticlePtr>::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<tPPtr,vector<LorentzRotation> >::const_iterator bit=_boosts.begin();bit!=_boosts.end();++bit) {
for(vector<LorentzRotation>::const_reverse_iterator rit=bit->second.rbegin();rit!=bit->second.rend();++rit) {
LorentzRotation rot = rit->inverse();
bit->first->transform(rot);
}
}
_boosts.clear();
for(map<tShowerTreePtr,vector<LorentzRotation> >::const_iterator bit=_treeBoosts.begin();bit!=_treeBoosts.end();++bit) {
for(vector<LorentzRotation>::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<ShowerParticlePtr>(p->parents()[0]);
if(parent) {
emitted=true;
reconstructSpaceLikeJet(parent);
}
// if branching reconstruct time-like child
if(p->children().size()==2)
child = dynamic_ptr_cast<ShowerParticlePtr>(p->children()[1]);
if(p->perturbative()==0 && child) {
dynamic_ptr_cast<ShowerParticlePtr>(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<ShowerProgenitorPtr> ShowerHardJets=decay->extractProgenitors();
for(unsigned int ix=0;ix<ShowerHardJets.size();++ix) {
_boosts[ShowerHardJets[ix]->progenitor()] = vector<LorentzRotation>();
}
for(map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator
tit = _currentTree->treelinks().begin();
tit != _currentTree->treelinks().end();++tit) {
_treeBoosts[tit->first] = vector<LorentzRotation>();
}
try {
bool radiated[2]={false,false};
// find the decaying particle and check if particles radiated
ShowerProgenitorPtr initial;
for(unsigned int ix=0;ix<ShowerHardJets.size();++ix) {
// only consider initial-state jets
if(ShowerHardJets[ix]->progenitor()->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;ix<ShowerHardJets.size();++ix) {
// only consider final-state jets
if(!ShowerHardJets[ix]->progenitor()->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<tPPtr,vector<LorentzRotation> >::const_iterator bit=_boosts.begin();bit!=_boosts.end();++bit) {
for(vector<LorentzRotation>::const_reverse_iterator rit=bit->second.rbegin();rit!=bit->second.rend();++rit) {
LorentzRotation rot = rit->inverse();
bit->first->transform(rot);
}
}
_boosts.clear();
for(map<tShowerTreePtr,vector<LorentzRotation> >::const_iterator bit=_treeBoosts.begin();bit!=_treeBoosts.end();++bit) {
for(vector<LorentzRotation>::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<tShowerTreePtr,pair<tShowerProgenitorPtr,
tShowerParticlePtr> >::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<tPPtr,vector<LorentzRotation> >::const_iterator bit=_boosts.begin();bit!=_boosts.end();++bit) {
for(vector<LorentzRotation>::const_reverse_iterator rit=bit->second.rbegin();rit!=bit->second.rend();++rit) {
LorentzRotation rot = rit->inverse();
bit->first->transform(rot);
}
}
_boosts.clear();
for(map<tShowerTreePtr,vector<LorentzRotation> >::const_iterator bit=_treeBoosts.begin();bit!=_treeBoosts.end();++bit) {
for(vector<LorentzRotation>::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<ShowerParticlePtr>(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<ShowerParticlePtr>(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<Energy2> pmag;
Energy total(Ejet);
for(unsigned int ix=0;ix<jetKinematics.size();++ix) {
pmag.push_back(jetKinematics[ix].p.vect().mag2());
total+=jetKinematics[ix].q.mass();
}
// return if no possible solution
if(total>mb) 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;iy<jetKinematics.size();++iy) {
if(jetKinematics[iy].parent==partner) continue;
ea = sqrt(sqr(d2)*pmag[iy]+jetKinematics[iy].q.mass2());
roots += ea;
ds += d2/ea*pmag[iy];
}
if(partner) {
ec = sqrt(sqr(d1)*pcmag + pt2 + ppartner[1].mass2());
roots += ec;
ds += d1/ec*pcmag;
}
d1 += (mb-roots)/ds;
d2 = d1 + pjn/pcn;
}
while(abs(mb-roots)>eps && 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<Lorentz5Momentum> pin;
vector<Lorentz5Momentum> pout;
// on-shell masses of the decay products
vector<Energy> mon;
Energy mbar(-GeV);
// the hard branchings of the particles
set<HardBranchingPtr>::iterator cit;
set<HardBranchingPtr> 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;ix<branch->children().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;ix<pout.size();++ix) pout[ix].boost(boostv);
if(initial->branchingParticle()->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;ix<pout.size();++ix) qtotal+=pout[ix];
Lorentz5Momentum qperp =
qisr-(qisr.vect()*qtotal.vect())/(qtotal.vect().mag2())*qtotal;
pvect +=qperp;
pvect /=k2;
pvect.setMass(mbar);
pvect.rescaleEnergy();
pvect.boost(-boostv);
(*cit)->pVector(pvect);
(*cit)->showerMomentum(pvect);
}
}
// For initial-state if needed
if(initial) {
tShowerParticlePtr newPartner=initial->branchingParticle()->partner();
if(newPartner) {
tHardBranchingPtr branch;
for( set<HardBranchingPtr>::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<HardBranchingPtr>::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<Lorentz5Momentum> pout,
vector<Energy> 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<Energy2> pmag;
for(unsigned int ix=0;ix<pout.size();++ix) {
pmag.push_back(pout[ix].vect().mag2());
}
// Newton-Raphson for the rescaling
vector<Energy> root(pout.size());
do {
// compute new energies
Energy sum(ZERO);
for(unsigned int ix=0;ix<pout.size();++ix) {
root[ix] = sqrt(pmag[ix]/sqr(lambda)+sqr(mon[ix]));
sum+=root[ix];
}
// if accuracy reached exit
if(abs(sum/roots-1.)<1e-10) break;
// use Newton-Raphson to compute new guess for lambda
Energy numer(ZERO),denom(ZERO);
for(unsigned int ix=0;ix<pout.size();++ix) {
numer +=root[ix];
denom +=pmag[ix]/root[ix];
}
numer-=roots;
double fact = 1.+sqr(lambda)*numer/denom;
if(fact<0.) fact=0.5;
lambda *=fact;
++ntry;
}
while(ntry<100);
}
if(std::isnan(lambda))
throw Exception() << "Rescaling factor is nan in KinematicsReconstructor::"
<< "inverseRescalingFactor "
<< Exception::eventerror;
return lambda;
}
bool KinematicsReconstructor::
deconstructGeneralSystem(HardTreePtr tree,
ShowerInteraction type) const {
// extract incoming and outgoing particles
ColourSingletShower in,out;
for(set<HardBranchingPtr>::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<HardBranchingPtr>::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<HardBranchingPtr>::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<ColourSingletShower>
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;ix<systems.size();++ix) {
if(systems[ix].type==II)
deconstructInitialInitialSystem(applyBoost,toRest,fromRest,tree,
systems[ix].jets,type);
}
if(type!=ShowerInteraction::QCD) {
combineFinalState(systems);
general=false;
}
}
// DIS and VBF type
else if(nnun==0&&nnii==0&&((nnif==1&&nnf>0&&nni==1)||
(nnif==2&& nni==0))) {
for(unsigned int ix=0;ix<systems.size();++ix) {
if(systems[ix].type==IF)
deconstructInitialFinalSystem(tree,systems[ix].jets,type);
}
}
// e+e- type
else if(nnun==0&&nnii==0&&nnif==0&&nnf>0&&nni==2) {
// only FS needed
// but need to boost to rest frame if QED ISR
Lorentz5Momentum ptotal;
for(unsigned int ix=0;ix<systems.size();++ix) {
if(systems[ix].type==I)
ptotal += systems[ix].jets[0]->branchingParticle()->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<systems.size();++ix) {
if(systems[ix].type==F)
deconstructFinalStateSystem(toRest,fromRest,tree,
systems[ix].jets,type);
}
// only at this point that we can be sure all the reference vectors
// are correct
for(set<HardBranchingPtr>::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<HardBranchingPtr>::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<HardBranchingPtr>::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<HardBranchingPtr>::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<HardBranchingPtr>::const_iterator it=tree->incoming().begin();
it!=tree->incoming().end();++it) {
(**it).setMomenta(LorentzRotation(),1.,Lorentz5Momentum(),false);
}
for(set<HardBranchingPtr>::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<ShowerProgenitorPtr> jets) const {
Lorentz5Momentum pin[2],pout[2],pbeam;
for(unsigned int ix=0;ix<jets.size();++ix) {
// final-state parton
if(jets[ix]->progenitor()->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;ix<jets.size();++ix) {
if(jets[ix]->progenitor()->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;ix<jets.size();++ix) {
if(jets[ix]->progenitor()->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<ShowerProgenitorPtr> jets) const {
bool added=false;
// add the intrinsic pt if needed
for(unsigned int ix=0;ix<jets.size();++ix) {
// only for initial-state particles which haven't radiated
if(jets[ix]->progenitor()->isFinalState()||
jets[ix]->hasEmitted()||
jets[ix]->reconstructed()==ShowerProgenitor::dontReconstruct) continue;
if(_intrinsic.find(jets[ix])==_intrinsic.end()) continue;
pair<Energy,double> 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<Energy2> 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<Energy2> 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()<ZERO ? -sqr(p.mass())+sqr(p.x())+sqr(p.y()) : sqr(p.mass())+sqr(p.x())+sqr(p.y()) ;
double ratio = mt2/(sqr(p.t())+sqr(q.t()));
if(abs(ratio)>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<ShowerProgenitorPtr> 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;ix<jets.size();++ix) {
radiated |=jets[ix]->hasEmitted();
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<ShowerProgenitorPtr>::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(num<ZERO) throw KinematicsReconstructionVeto();
double k = sqrt( num / (4.*s*jit->p.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<ShowerProgenitorPtr> 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;ix<jets.size();++ix) {
if(jets[ix]->reconstructed()==ShowerProgenitor::notReconstructed)
jets[ix]->reconstructed(ShowerProgenitor::done);
}
return;
}
// initial state shuffling
applyBoost=false;
vector<Lorentz5Momentum> p, pq, p_in;
vector<Energy> pts;
for(unsigned int ix=0;ix<jets.size();++ix) {
// add momentum to vector
p_in.push_back(jets[ix]->progenitor()->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;ix<jets.size();++ix) {
p.push_back(jets[ix]->progenitor()->momentum());
}
double x1 = p_in[0].z()/pq[0].z();
double x2 = p_in[1].z()/pq[1].z();
vector<double> 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;ix<jets.size();++ix) {
tPPtr toBoost = jets[ix]->progenitor();
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()<ZERO||newcmf.e()<ZERO) throw KinematicsReconstructionVeto();
findInitialBoost(pcm,newcmf,toRest,fromRest);
}
void KinematicsReconstructor::
deconstructInitialInitialSystem(bool & applyBoost,
LorentzRotation & toRest,
LorentzRotation & fromRest,
HardTreePtr tree,
vector<HardBranchingPtr> jets,
ShowerInteraction) const {
assert(jets.size()==2);
// put beam with +z first
if(jets[0]->beam()->momentum().z()<ZERO) swap(jets[0],jets[1]);
// get the momenta of the particles
vector<Lorentz5Momentum> pin,pq;
for(unsigned int ix=0;ix<jets.size();++ix) {
pin.push_back(jets[ix]->branchingParticle()->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<double> boost = inverseInitialStateRescaling(x[0],x[1],pcm,pin,pq);
set<HardBranchingPtr>::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<HardBranchingPtr> 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<Lorentz5Momentum> ptemp;
set<HardBranchingPtr>::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<tcQTildeShowerHandlerPtr>(ShowerHandler::currentHandler())->partnerFinder()
->setInitialEvolutionScales(particles,false,type,false);
// calculate the reference vectors
unsigned int iloc(0);
set<HardBranchingPtr>::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<HardBranchingPtr>::iterator cit;
vector<Lorentz5Momentum> pout;
vector<Energy> 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;ix<pout.size();++ix) {
pout[ix].transform(trans);
}
// rescaling factor
double lambda=inverseRescalingFactor(pout,mon,pin.mass());
if (lambda< 1.e-10) throw KinematicsReconstructionVeto();
// now calculate the p reference vectors
for(unsigned int ix=0;ix<jets.size();++ix) {
Lorentz5Momentum pvect = jets[ix]->branchingParticle()->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<Lorentz5Momentum> ptemp;
set<HardBranchingPtr>::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<tcQTildeShowerHandlerPtr>(ShowerHandler::currentHandler())->partnerFinder()
->setInitialEvolutionScales(particles,false,type,false);
// calculate the reference vectors
unsigned int iloc(0);
set<HardBranchingPtr>::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<ShowerParticlePtr>(p->children()[1]))
deepTransform(p->children()[1],bv);
}
}
namespace {
bool sortJets(ShowerProgenitorPtr j1, ShowerProgenitorPtr j2) {
return j1->highestpT()>j2->highestpT();
}
}
void KinematicsReconstructor::
reconstructGeneralSystem(vector<ShowerProgenitorPtr> & ShowerHardJets) const {
// find initial- and final-state systems
ColourSingletSystem in,out;
for(unsigned int ix=0;ix<ShowerHardJets.size();++ix) {
if(ShowerHardJets[ix]->progenitor()->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<ShowerProgenitorPtr> & ShowerHardJets) const {
static const Energy2 minQ2 = 1e-4*GeV2;
map<ShowerProgenitorPtr,bool> used;
for(unsigned int ix=0;ix<ShowerHardJets.size();++ix) {
used[ShowerHardJets[ix]] = false;
} // first to the final-state reconstruction of any systems which need it
set<ShowerProgenitorPtr> outgoing;
// first find any particles with final state partners
for(unsigned int ix=0;ix<ShowerHardJets.size();++ix) {
if(ShowerHardJets[ix]->progenitor()->isFinalState()&&
ShowerHardJets[ix]->progenitor()->partner()&&
ShowerHardJets[ix]->progenitor()->partner()->isFinalState()) outgoing.insert(ShowerHardJets[ix]);
}
// then find the colour partners
if(!outgoing.empty()) {
set<ShowerProgenitorPtr> partners;
for(set<ShowerProgenitorPtr>::const_iterator it=outgoing.begin();it!=outgoing.end();++it) {
for(unsigned int ix=0;ix<ShowerHardJets.size();++ix) {
if((**it).progenitor()->partner()==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<ShowerProgenitorPtr> outgoingJets(outgoing.begin(),outgoing.end());
reconstructFinalStateSystem(false,toRest,fromRest,outgoingJets);
}
// Now do any initial-final systems which are needed
vector<ColourSingletSystem> 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;ix<ShowerHardJets.size();++ix) {
if(!ShowerHardJets[ix]->progenitor()->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;is<IFSystems.size();++is) {
for(unsigned int ix=0;ix<ShowerHardJets.size();++ix) {
if(IFSystems[is].jets[0]->progenitor()->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;ix<IFSystems[is].jets.size();++ix) {
if(IFSystems[is].jets[ix]->progenitor()->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;ix<IFSystems[is].jets.size();++ix) {
used[IFSystems[is].jets[ix]] = true;
}
}
}
}
// now we finally need to handle the initial state system
ColourSingletSystem in,out;
for(unsigned int ix=0;ix<ShowerHardJets.size();++ix) {
if(ShowerHardJets[ix]->progenitor()->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<in.jets.size();++ix) {
if(!used[in.jets[ix]]) {
doRecon = true;
break;
}
}
LorentzRotation toRest,fromRest;
bool applyBoost(false);
if(doRecon) {
reconstructInitialInitialSystem(applyBoost,toRest,fromRest,in.jets);
}
// reconstruct the final-state systems
if(!doRecon) {
for(unsigned int ix=0;ix<out.jets.size();++ix) {
if(!used[out.jets[ix]]) {
doRecon = true;
break;
}
}
}
if(doRecon) {
reconstructFinalStateSystem(applyBoost,toRest,fromRest,out.jets);
}
}
void KinematicsReconstructor::
reconstructColourPartner(vector<ShowerProgenitorPtr> & 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<ShowerParticlePtr,ShowerProgenitorPtr> progenitorMap;
for(unsigned int ix=0;ix<ShowerHardJets.size();++ix) {
progenitorMap[ShowerHardJets[ix]->progenitor()] = ShowerHardJets[ix];
}
// check that the IF systems can be reconstructed
bool canReconstruct = true;
for(unsigned int ix=0;ix<ShowerHardJets.size();++ix) {
tShowerParticlePtr progenitor = ShowerHardJets[ix]->progenitor();
tShowerParticlePtr partner = progenitor->partner();
if(!partner) continue;
else if((progenitor->isFinalState() &&
!partner->isFinalState()) ||
(!progenitor->isFinalState() &&
partner->isFinalState()) ) {
vector<ShowerProgenitorPtr> jets(2);
jets[0] = ShowerHardJets[ix];
jets[1] = progenitorMap[partner];
Lorentz5Momentum psum;
for(unsigned int iy=0;iy<jets.size();++iy) {
if(jets[iy]->progenitor()->isFinalState())
psum += jets[iy]->progenitor()->momentum();
else
psum -= jets[iy]->progenitor()->momentum();
}
if(-psum.m2()<minQ2) {
canReconstruct = false;
break;
}
}
}
if(!canReconstruct) {
reconstructGeneralSystem(ShowerHardJets);
return;
}
map<ShowerProgenitorPtr,bool> used;
for(unsigned int ix=0;ix<ShowerHardJets.size();++ix) {
used[ShowerHardJets[ix]] = false;
}
for(unsigned int ix=0;ix<ShowerHardJets.size();++ix) {
// skip jets which have already been handled
if(ShowerHardJets[ix]->reconstructed()==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<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::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<ShowerProgenitorPtr> 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<ShowerProgenitorPtr> 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;iy<jets.size();++iy) {
if(jets[iy]->progenitor()->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;iy<ShowerHardJets.size();++iy) {
if(ShowerHardJets[iy]->progenitor()->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;iy<out.jets.size();++iy) {
if(out.jets[iy]->reconstructed()==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;iy<out.jets.size();++iy) {
deepTransform(out.jets[iy]->progenitor(),finalBoosts);
}
for(unsigned int iy=0;iy<out.jets.size();++iy) {
if(out.jets[iy]->reconstructed()==ShowerProgenitor::dontReconstruct)
out.jets[iy]->reconstructed(ShowerProgenitor::notReconstructed);
}
}
}
}
bool KinematicsReconstructor::
inverseDecayRescalingFactor(vector<Lorentz5Momentum> pout,
vector<Energy> mon,Energy roots,
Lorentz5Momentum ppartner, Energy mbar,
double & k1, double & k2) const {
ThreeVector<Energy> qtotal;
vector<Energy2> pmag;
for(unsigned int ix=0;ix<pout.size();++ix) {
pmag.push_back(pout[ix].vect().mag2());
qtotal+=pout[ix].vect();
}
Energy2 dot1 = qtotal*ppartner.vect();
Energy2 qmag2=qtotal.mag2();
double a = -dot1/qmag2;
static const Energy eps=1e-10*GeV;
unsigned int itry(0);
Energy numer(ZERO),denom(ZERO);
k1=1.;
do {
++itry;
numer=denom=0.*GeV;
double k12=sqr(k1);
for(unsigned int ix=0;ix<pout.size();++ix) {
Energy en = sqrt(pmag[ix]/k12+sqr(mon[ix]));
numer += en;
denom += pmag[ix]/en;
}
Energy en = sqrt(qmag2/k12+sqr(mbar));
numer += en-roots;
denom += qmag2/en;
k1 += numer/denom*k12*k1;
if(abs(k1)>1e10) return false;
}
while (abs(numer)>eps&&itry<100);
k1 = abs(k1);
k2 = a*k1;
return itry<100;
}
void KinematicsReconstructor::
deconstructInitialFinalSystem(HardTreePtr tree,vector<HardBranchingPtr> jets,
ShowerInteraction type) const {
HardBranchingPtr incoming;
Lorentz5Momentum pin[2],pout[2],pbeam;
HardBranchingPtr initial;
Energy mc(ZERO);
for(unsigned int ix=0;ix<jets.size();++ix) {
// final-state parton
if(jets[ix]->status()==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;ix<jets.size();++ix) {
if(jets[ix]->status()==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<Lorentz5Momentum> ptemp;
set<HardBranchingPtr>::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<tcQTildeShowerHandlerPtr>(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<HardBranchingPtr>::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<HardBranchingPtr>::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<HardBranchingPtr>::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<ShowerParticlePtr>(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<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::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(num<ZERO) throw KinematicsReconstructionVeto();
double k = sqrt( num / (4.*s*orderedJets.begin()->p.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;ix<branch->children().size();++ix) {
if(branch->branchingParticle()->id()==
branch->children()[ix]->branchingParticle()->id())
return findMass(branch->children()[ix]);
}
}
return branch->branchingParticle()->dataPtr()->mass();
}
vector<double>
KinematicsReconstructor::inverseInitialStateRescaling(double & x1, double & x2,
const Lorentz5Momentum & pold,
const vector<Lorentz5Momentum> & p,
const vector<Lorentz5Momentum> & 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<double> 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<double>
KinematicsReconstructor::initialStateRescaling(double x1, double x2,
const Lorentz5Momentum & pold,
const vector<Lorentz5Momentum> & p,
const vector<Lorentz5Momentum> & pq,
const vector<Energy>& 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()<ZERO ? exp(-2.*pold.rapidity()) : exp(2.*pold.rapidity());
if(rad <= 0.) throw KinematicsReconstructionVeto();
k1 = sqrt(rad);
k2 = kp/k1;
}
// conserve longitudinal momentum
else if(_initialStateReconOption==1) {
double a2 = (a[0]+a[1]/kp);
double b2 = -x2+x1;
double c2 = -(b[1]*kp+b[0]);
if(abs(b2)>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<double> 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<ShowerProgenitorPtr> & ShowerHardJets,
ShowerInteraction type) const {
// identify and catagorize the colour singlet systems
unsigned int nnun(0),nnii(0),nnif(0),nnf(0),nni(0);
vector<ColourSingletSystem>
systems(identifySystems(set<ShowerProgenitorPtr>(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;ix<systems.size();++ix) {
if(systems[ix].type==II)
reconstructInitialInitialSystem(applyBoost,toRest,fromRest,
systems[ix].jets);
}
if(type!=ShowerInteraction::QCD) {
combineFinalState(systems);
general=false;
}
}
// DIS and VBF type
else if(nnun==0&&nnii==0&&((nnif==1&&nnf>0&&nni==1)||
(nnif==2&& nni==0))) {
// check these systems can be reconstructed
for(unsigned int ix=0;ix<systems.size();++ix) {
// compute q^2
if(systems[ix].type!=IF) continue;
Lorentz5Momentum q;
for(unsigned int iy=0;iy<systems[ix].jets.size();++iy) {
if(systems[ix].jets[iy]->progenitor()->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;ix<systems.size();++ix) {
if(systems[ix].type==IF)
reconstructInitialFinalSystem(systems[ix].jets);
}
}
}
// e+e- type
else if(nnun==0&&nnii==0&&nnif==0&&nnf>0&&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<systems.size();++ix) {
if(systems[ix].type==F)
reconstructFinalStateSystem(applyBoost,toRest,fromRest,
systems[ix].jets);
}
}
else {
reconstructGeneralSystem(ShowerHardJets);
}
}
void KinematicsReconstructor::findInitialBoost(const Lorentz5Momentum & pold,
const Lorentz5Momentum & pnew,
LorentzRotation & toRest,
LorentzRotation & fromRest) const {
// do one boost
if(_initialBoost==0) {
toRest = LorentzRotation(pold.findBoostToCM());
fromRest = LorentzRotation(pnew.boostVector());
}
else if(_initialBoost==1) {
// boost to rest frame
// first transverse
toRest = Boost(-pold.x()/pold.t(),-pold.y()/pold.t(),0.);
// then longitudinal
double beta = pold.z()/sqrt(pold.m2()+sqr(pold.z()));
toRest.boost((Boost(0.,0.,-beta)));
// boost from rest frame
// first apply longitudinal boost
beta = pnew.z()/sqrt(pnew.m2()+sqr(pnew.z()));
fromRest=LorentzRotation(Boost(0.,0.,beta));
// then transverse one
fromRest.boost(Boost(pnew.x()/pnew.t(),
pnew.y()/pnew.t(),0.));
}
else
assert(false);
}
diff --git a/Shower/QTilde/Kinematics/KinematicsReconstructor.h b/Shower/QTilde/Kinematics/KinematicsReconstructor.h
--- a/Shower/QTilde/Kinematics/KinematicsReconstructor.h
+++ b/Shower/QTilde/Kinematics/KinematicsReconstructor.h
@@ -1,656 +1,655 @@
// -*- C++ -*-
//
// KinematicsReconstructor.h 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.
//
#ifndef HERWIG_KinematicsReconstructor_H
#define HERWIG_KinematicsReconstructor_H
//
// This is the declaration of the KinematicsReconstructor class.
//
#include "ThePEG/Interface/Interfaced.h"
#include "Herwig/Shower/QTilde/Base/ShowerParticle.h"
#include "Herwig/Shower/QTilde/Base/ShowerProgenitor.h"
#include "Herwig/Shower/QTilde/Base/ShowerTree.h"
#include "Herwig/Shower/QTilde/Base/HardTree.h"
#include "KinematicsReconstructor.fh"
#include <cassert>
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<JetKinStruct> JetKinVect;
/**
* Enum to identify types of colour singlet systems
*/
enum SystemType { UNDEFINED=-1, II, IF, F ,I };
/**
* Struct to store colour singlets
*/
template<typename Value> struct ColourSinglet {
typedef vector<ColourSinglet<Value> > 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<Value> jets;
};
/**
* Struct to store a colour singlet system of particles
*/
typedef ColourSinglet<ShowerProgenitorPtr> ColourSingletSystem;
/**
* Struct to store a colour singlet shower
*/
typedef ColourSinglet<HardBranchingPtr> 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<tShowerProgenitorPtr,
pair<Energy,double> > & 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<ShowerProgenitorPtr>) const;
/**
* Perform the reconstruction of a system with only final-state
* particles
*/
void reconstructFinalStateSystem(bool applyBoost,
const LorentzRotation & toRest,
const LorentzRotation & fromRest,
vector<ShowerProgenitorPtr>) const;
/**
* Reconstruction of a general coloured system
*/
void reconstructGeneralSystem(vector<ShowerProgenitorPtr> & ShowerHardJets) const;
/**
* Reconstruction of a general coloured system doing
* final-final, then initial-final and then initial-initial
*/
void reconstructFinalFirst(vector<ShowerProgenitorPtr> & ShowerHardJets) const;
/**
* Reconstruction of a general coloured system doing
* colour parners
*/
void reconstructColourPartner(vector<ShowerProgenitorPtr> & ShowerHardJets) const;
/**
* Reconstruction based on colour singlet systems
*/
void reconstructColourSinglets(vector<ShowerProgenitorPtr> & ShowerHardJets,
ShowerInteraction type) const;
/**
* Perform the reconstruction of a system with only final-state
* particles
*/
void reconstructInitialInitialSystem(bool & applyBoost,
LorentzRotation & toRest,
LorentzRotation & fromRest,
vector<ShowerProgenitorPtr>) 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<HardBranchingPtr>,
ShowerInteraction) const;
/**
* Perform the inverse reconstruction of a system with only initial-state
* particles
*/
void deconstructInitialInitialSystem(bool & applyBoost,
LorentzRotation & toRest,
LorentzRotation & fromRest,
HardTreePtr,
vector<HardBranchingPtr>,
ShowerInteraction ) const;
/**
* Perform the inverse reconstruction of a system with only initial-state
* particles
*/
void deconstructInitialFinalSystem(HardTreePtr,
vector<HardBranchingPtr>,
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<Lorentz5Momentum> pout,
vector<Energy> 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<Lorentz5Momentum> pout,
vector<Energy> 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<typename Value> void findPartners(Value branch,set<Value> & done,
const set<Value> & branchings,
vector<Value> & jets) const;
/**
* Add the intrinsic \f$p_T\f$ to the system if needed
*/
bool addIntrinsicPt(vector<ShowerProgenitorPtr>) 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<double> initialStateRescaling(double x1, double x2,
const Lorentz5Momentum & pold,
const vector<Lorentz5Momentum> & p,
const vector<Lorentz5Momentum> & pq,
const vector<Energy>& highespts) const;
/**
* Calculate the inverse of the initial-state rescaling factor
*/
vector<double> inverseInitialStateRescaling(double & x1, double & x2,
const Lorentz5Momentum & pold,
const vector<Lorentz5Momentum> & p,
const vector<Lorentz5Momentum> & pq) const;
/**
* Find the colour singlet systems
*/
template<typename Value >
typename ColourSinglet<Value>::VecType identifySystems(set<Value> jets,
unsigned int & nnun,unsigned int & nnii,
unsigned int & nnif,unsigned int & nnf,
unsigned int & nni) const;
/**
* Combine final-state colour systems
*/
template<typename Value>
void combineFinalState(vector<ColourSinglet<Value> > & 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<tShowerProgenitorPtr,pair<Energy,double> > _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<cPDPtr> _noRescale;
/**
* Storage of the boosts applied to enable resetting after failure
*/
mutable map<tPPtr,vector<LorentzRotation> > _boosts;
/**
* Storage of the boosts applied to enable resetting after failure
*/
mutable map<tShowerTreePtr,vector<LorentzRotation> > _treeBoosts;
};
}
-#include "KinematicsReconstructor.tcc"
#endif /* HERWIG_KinematicsReconstructor_H */

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 5:57 PM (9 h, 36 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023718
Default Alt Text
(244 KB)

Event Timeline