diff --git a/.gitignore b/.gitignore index 0b8abd75..165ebe0b 100644 --- a/.gitignore +++ b/.gitignore @@ -2,8 +2,11 @@ *.pyc __init__.py +MELA/data/*.tar.gz* +MELA/data/*/*.tar.gz* MELA/data/*/libmcfm*.so MELA/data/*/libjhugen*.so +MELA/data/*/libcollier*.so MELA/src/*.o MELA/src/*.so @@ -13,6 +16,18 @@ MELA/interface/*.o MELA/interface/*.so MELA/interface/*.d +MELA/COLLIER/*.o +MELA/COLLIER/*.mod +MELA/COLLIER/*.so +MELA/COLLIER/*/* +MELA/COLLIER/*.f +MELA/COLLIER/*.F +MELA/COLLIER/*.F90 +MELA/COLLIER/*.c +MELA/COLLIER/*.cc +MELA/COLLIER/*.h +MELA/COLLIER/*.hh + MELA/fortran/*.o MELA/fortran/*.mod MELA/fortran/*.so @@ -29,7 +44,8 @@ MELA/test/*.out MELA/test/*.ref MELA/test/reference/*.out MELA/test/*.o -MELA/test/*_c.so -MELA/test/*_c.d +MELA/test/*.a +MELA/test/*.so +MELA/test/*.d MELA/test/*.pcm diff --git a/MELA/BuildFile.xml b/MELA/BuildFile.xml index 58470fa7..6f32c526 100755 --- a/MELA/BuildFile.xml +++ b/MELA/BuildFile.xml @@ -4,7 +4,7 @@ - + @@ -22,3 +22,8 @@ + + + + + diff --git a/MELA/COLLIER/makefile b/MELA/COLLIER/makefile new file mode 100644 index 00000000..2010a3a1 --- /dev/null +++ b/MELA/COLLIER/makefile @@ -0,0 +1,118 @@ +MAINDIR = . +AUXDIR = $(MAINDIR)/Aux +COLIDIR = $(MAINDIR)/COLI +TENSORSDIR = $(MAINDIR)/tensors +DDLIBDIR = $(MAINDIR)/DDlib +#MODULES=$(MAINDIR)/modules +#OBJECTS=$(MAINDIR)/objects +MODULES=$(MAINDIR) +OBJECTS=$(MAINDIR) +ROOFITINCLUDE = +RM = /bin/rm +INCLUDE = -I$(MAINDIR) $(ROOFITINCLUDE) +Comp = gfort + +ifeq ($(Comp),ifort) + fcomp = ifort -132 -O2-lifcore -fPIC +endif +ifeq ($(Comp),gfort) + fcomp = gfortran -ffixed-line-length-132 -fno-default-integer-8 -fPIC -O2 -funroll-loops -Wtabs +endif + +ifeq ($(Comp),gfort) + ccomp = gcc +endif +ifeq ($(Comp),ifort) + ccomp = icc +endif + +ifeq ($(Comp),gfort) + clib = -lm -lgfortran +endif +ifeq ($(Comp),ifort) + clib = -lm -lirc +endif + + +LIB = libcollier.so + +Deps = $(MAINDIR)/collier_global.F90 \ + $(AUXDIR)/Combinatorics.F90 \ + $(AUXDIR)/master.F90 \ + $(COLIDIR)/coli_aux2.F90 \ + $(COLIDIR)/coli_stat.F90 \ + $(MAINDIR)/collier_aux.F90 \ + $(AUXDIR)/cache.F90 \ + $(TENSORSDIR)/InitTensors.F90 \ + $(MAINDIR)/collier_init.F90 \ + $(COLIDIR)/reductionAB.F90 \ + $(COLIDIR)/reductionC.F90 \ + $(COLIDIR)/reductionD.F90 \ + $(COLIDIR)/reductionEFG.F90 \ + $(COLIDIR)/reductionTN.F90 \ + $(MAINDIR)/collier_coefs.F90 \ + $(TENSORSDIR)/BuildTensors.F90 \ + $(TENSORSDIR)/TensorReduction.F90 \ + $(MAINDIR)/collier_tensors.F90 \ + $(DDLIBDIR)/DD_global.F90 \ + $(MAINDIR)/COLLIER.F90 \ + $(COLIDIR)/coli_b0.F \ + $(COLIDIR)/coli_d0.F \ + $(COLIDIR)/coli_d0reg.F \ + $(COLIDIR)/coli_aux.F \ + $(COLIDIR)/coli_c0.F \ + $(DDLIBDIR)/DD_aux.F \ + $(DDLIBDIR)/DD_5pt.F \ + $(DDLIBDIR)/DD_6pt.F \ + $(DDLIBDIR)/DD_3pt.F \ + $(DDLIBDIR)/DD_to_COLLIER.F \ + $(DDLIBDIR)/DD_2pt.F + +Objs = collier_global.o \ + Aux_Combinatorics.o \ + Aux_master.o \ + COLI_coli_aux2.o \ + COLI_coli_stat.o \ + collier_aux.o \ + Aux_cache.o \ + tensors_InitTensors.o \ + collier_init.o \ + COLI_reductionAB.o \ + COLI_reductionC.o \ + COLI_reductionD.o \ + COLI_reductionEFG.o \ + COLI_reductionTN.o \ + collier_coefs.o \ + tensors_BuildTensors.o \ + tensors_TensorReduction.o \ + collier_tensors.o \ + DDlib_DD_global.o \ + COLLIER.o \ + COLI_coli_b0.o \ + COLI_coli_d0.o \ + COLI_coli_d0reg.o \ + COLI_coli_aux.o \ + COLI_coli_c0.o \ + DDlib_DD_aux.o \ + DDlib_DD_5pt.o \ + DDlib_DD_6pt.o \ + DDlib_DD_3pt.o \ + DDlib_DD_to_COLLIER.o \ + DDlib_DD_2pt.o + + +${Objs}: $(Deps) + @echo " " + @echo " Compiling COLLIER dependencies" + $(fcomp) -c $(Deps) + @echo " " + @echo " Compiling COLLIER library" + g++ -Wl,-soname,$(LIB) -shared -o $(LIB) *.o + +clean: + @echo " deleting object files" + rm -f *.so $(OBJECTS)/*.o $(MODULES)/*.mod + + +# supresses command calls +.SILENT: diff --git a/MELA/COLLIER/setup.sh b/MELA/COLLIER/setup.sh new file mode 100644 index 00000000..575e9c3d --- /dev/null +++ b/MELA/COLLIER/setup.sh @@ -0,0 +1,49 @@ +#!/bin/bash + +{ + +scriptdir=$(dirname $0) +curdir=$(pwd) + +cd $scriptdir + +pkgname="collier-1.2" +pkgdir="COLLIER-1.2" +tarname=$pkgname".tar.gz" +tarweb="https://www.hepforge.org/archive/collier/"$tarname +libname="libcollier.so" +tmpdir="colliertmp" + +if [[ $# > 0 ]] && [[ "$1" == *"clean"* ]];then + + rm -f *.so + rm -f *.o + rm -f *.mod + rm -f *.f + rm -f *.F + rm -f *.F90 + for f in $(ls ./);do + if [ -d $f ];then + rm -rf $f + fi + done + + rm -f "../data/"$SCRAM_ARCH"/"$libname + +else + + wget $tarweb + mkdir $tmpdir + tar -xvzf $tarname -C $tmpdir + rm $tarname + mv $tmpdir"/"$pkgdir"/src/"* ./ + rm -rf $tmpdir + + make + mv $libname "../data/"$SCRAM_ARCH"/"$libname + +fi + +cd $curdir + +} \ No newline at end of file diff --git a/MELA/fortran/makefile b/MELA/fortran/makefile index eb93a0ea..8a00d658 100755 --- a/MELA/fortran/makefile +++ b/MELA/fortran/makefile @@ -4,14 +4,14 @@ all: false else -CMSROOT = ./ -#MODULES=$(CMSROOT)modules/ -#OBJECTS=$(CMSROOT)objects/ -MODULES=$(CMSROOT) -OBJECTS=$(CMSROOT) +MAINDIR = ./ +#MODULES=$(MAINDIR)modules/ +#OBJECTS=$(MAINDIR)objects/ +MODULES=$(MAINDIR) +OBJECTS=$(MAINDIR) ROOFITINCLUDE = RM = /bin/rm -INCLUDE = -I$(CMSROOT) $(ROOFITINCLUDE) +INCLUDE = -I$(MAINDIR) $(ROOFITINCLUDE) Comp = gfort # Linking LHAPDF @@ -28,12 +28,26 @@ else LHAPDFflags = -DuseLHAPDF=0 endif +# Linking the Collier library +UseCOLLIER=Yes +# directory which contains libLHAPDF.a, libLHAPDF.la, libLHAPDF.so +MyCOLLIERDir=../data/${SCRAM_ARCH}/ +MyCOLLIERInc=../COLLIER/ +# remember to export +# LD_LIBRARY_PATH=/.../LHAPDF-x.y.z/lib/:${LD_LIBRARY_PATH} +# LHAPDF_DATA_PATH=/.../LHAPDF-x.y.z/share/LHAPDF/:${LHAPDF_DATA_PATH} +ifeq ($(UseCOLLIER),Yes) + COLLIERflags = -L$(MyCOLLIERDir) -lcollier -DuseCollier=1 -I$(MyCOLLIERInc) +else + COLLIERflags = -DuseCollier=0 +endif + ifeq ($(Comp),ifort) - fcomp = ifort -fpp -O2 -vec-report0 -Dcompiler=1 -lifcore -fPIC + fcomp = ifort -fpp -O2 -vec-report0 -Dcompiler=1 -lifcore $(LHAPDFflags) $(COLLIERflags) -fPIC endif ifeq ($(Comp),gfort) - fcomp = gfortran -O0 -ffree-line-length-none -Dcompiler=2 $(LHAPDFflags) -fno-automatic -fno-f2c -fPIC -g + fcomp = gfortran -O0 -ffree-line-length-none -Dcompiler=2 $(LHAPDFflags) $(COLLIERflags) -fno-automatic -fno-f2c -fPIC -g endif ifeq ($(Comp),gfort) diff --git a/MELA/fortran/mod_JHUGen.F90 b/MELA/fortran/mod_JHUGen.F90 index 478e7c44..d239b04a 100755 --- a/MELA/fortran/mod_JHUGen.F90 +++ b/MELA/fortran/mod_JHUGen.F90 @@ -26,6 +26,7 @@ SUBROUTINE InitFirstTime(pdftable,pdfstrlength,pdfmember) includeInterference=.true. includeGammaStar=.true. + MPhotonCutoff=4d0*GeV WidthScheme=0 PDFSet=3 ! 1: CTEQ6L1 2: MRSW with best fit, 2xx: MSTW with eigenvector set xx=01..40 LHAPDFString = pdftable @@ -81,6 +82,13 @@ SUBROUTINE InitFirstTime(pdftable,pdfstrlength,pdfmember) print *,"Collider not implemented." stop ENDIF + +#if useCollier==1 + Collier_maxNLoopProps = -1 + Collier_maxRank = -1 +#endif + call InitCOLLIER(4,3) ! Arguments for ggZH + return END SUBROUTINE @@ -175,7 +183,7 @@ SUBROUTINE InitPDFs() implicit none DOUBLE PRECISION alphasPDF - call InitPDFset(trim(LHAPDFString,lenLHAPDFString)) ! Let LHAPDF handle everything + call InitPDFSetByName(trim(LHAPDFString,lenLHAPDFString)) ! Let LHAPDF handle everything call InitPDF(LHAPDFMember) alphas_mz=alphasPDF(zmass_pdf) diff --git a/MELA/fortran/mod_Misc.F90 b/MELA/fortran/mod_Misc.F90 index 92fc65af..d30a69bd 100755 --- a/MELA/fortran/mod_Misc.F90 +++ b/MELA/fortran/mod_Misc.F90 @@ -14,6 +14,7 @@ MODULE ModMisc MODULE PROCEDURE MinkowskyProduct MODULE PROCEDURE MinkowskyProductC MODULE PROCEDURE MinkowskyProductRC + MODULE PROCEDURE MinkowskyProductCR END INTERFACE OPERATOR (.dot.) INTERFACE OPERATOR (.cross.) @@ -72,6 +73,18 @@ FUNCTION MinkowskyProductRC(p1,p2) END FUNCTION MinkowskyProductRC +FUNCTION MinkowskyProductCR(p1,p2) +implicit none +real(8), intent(in) :: p2(1:4) +complex(8), intent(in) :: p1(1:4) +complex(8) :: MinkowskyProductCR + + MinkowskyProductCR = p1(1)*p2(1) & + - p1(2)*p2(2) & + - p1(3)*p2(3) & + - p1(4)*p2(4) +END FUNCTION MinkowskyProductCR + double complex function et1(e1,e2,e3,e4) implicit none complex(8), intent(in) :: e1(4), e2(4), e3(4), e4(4) @@ -399,8 +412,6 @@ FUNCTION LeviCiv(e1,e2,e3,e4) END FUNCTION - - FUNCTION pol_gluon_incoming(p,hel) implicit none complex(8) :: pol_gluon_incoming(1:4) @@ -420,8 +431,6 @@ FUNCTION pol_gluon_incoming(p,hel) cphi= p(2)/pv/st sphi= p(3)/pv/st endif - - ! -- distinguish between positive and negative energies ! if ( p(1) .gt. 0.0_dp) then ! hel=hel @@ -437,15 +446,9 @@ FUNCTION pol_gluon_incoming(p,hel) pol_gluon_incoming(2) = ct*cphi - (0d0,1d0)*hel*sphi pol_gluon_incoming(3) = ct*sphi + (0d0,1d0)*hel*cphi pol_gluon_incoming(4) = -st - - pol_gluon_incoming(2:4) = pol_gluon_incoming(2:4)/dsqrt(2.0d0) - -END FUNCTION - - - - + pol_gluon_incoming(2:4) = pol_gluon_incoming(2:4)/dsqrt(2.0d0) +END FUNCTION FUNCTION pol_Zff_outgoing(pf,pfbar,hel) ! does not contain the division by 1/q^2 as in pol_dk2mom implicit none @@ -453,8 +456,6 @@ FUNCTION pol_Zff_outgoing(pf,pfbar,hel) ! does not contain the division by 1/q^ integer :: hel real(8):: pf(1:4),pfbar(1:4) complex(8) :: Ub(1:4),V(1:4) - - Ub(1:4) = ubar_spinor(pf,hel) V(1:4) = v_spinor(pfbar,-hel) @@ -463,9 +464,8 @@ FUNCTION pol_Zff_outgoing(pf,pfbar,hel) ! does not contain the division by 1/q^ pol_Zff_outgoing(2)=-(-Ub(1)*V(4)+V(1)*Ub(4)-Ub(2)*V(3)+V(2)*Ub(3)) pol_Zff_outgoing(3)=-(0d0,1d0)*(Ub(1)*V(4)+V(1)*Ub(4)-Ub(2)*V(3)-V(2)*Ub(3)) pol_Zff_outgoing(4)=-(Ub(2)*V(4)-V(2)*Ub(4)-Ub(1)*V(3)+V(1)*Ub(3)) - RETURN -END FUNCTION +END FUNCTION ! ubar spinor, massless @@ -509,8 +509,6 @@ FUNCTION ubar_spinor(p,hel) RETURN END FUNCTION - - ! -- v0 spinor, massless FUNCTION v_spinor(p,hel) implicit none @@ -552,6 +550,165 @@ FUNCTION v_spinor(p,hel) RETURN END FUNCTION +function Chir_Weyl(sign,sp) ! Chir = sp.omega_sign = omega_sign.sp +implicit none +integer :: sign +double complex :: sp(1:4) +double complex :: Chir_Weyl(1:4) + if(sign.eq.+1) then !omega_+ + Chir_Weyl(1) = sp(1) + Chir_Weyl(2) = sp(2) + Chir_Weyl(3) = 0d0 + Chir_Weyl(4) = 0d0 + else !omega_- + Chir_Weyl(1) = 0d0 + Chir_Weyl(2) = 0d0 + Chir_Weyl(3) = sp(3) + Chir_Weyl(4) = sp(4) + endif + return +end function + +FUNCTION vbqq_Weyl(sp1,sp2) +implicit none +complex(8), intent(in) :: sp1(:), sp2(:) +integer, parameter :: Dv=4 +integer :: i +complex(8) :: vbqq_Weyl(Dv) +complex(8) :: rr, va(Dv),sp1a(4) + + va=(0d0,0d0) + vbqq_Weyl=(0d0,0d0) + + do i=1,Dv + if (i.eq.1) then + va(1)=(1d0,0d0) + else + va(i)=(-1d0,0d0) + endif + sp1a=spb2_Weyl(sp1,va) + + + rr=sum(sp1a(1:4)*sp2(1:4)) + if (i.eq.1) then + vbqq_Weyl = vbqq_Weyl + rr*va + else + vbqq_Weyl = vbqq_Weyl - rr*va + endif + va(i)=(0d0,0d0) + enddo + +END FUNCTION + +function spb2_Weyl(sp,v) +implicit none +complex(8), intent(in) :: sp(:),v(:) +complex(8) :: spb2_Weyl(4) +complex(8) :: x0(4,4),xx(4,4),xy(4,4) +complex(8) :: xz(4,4),x5(4,4) +complex(8) :: y1,y2,y3,y4,bp,bm,cp,cm +integer :: i,i1,i2,i3,Dv,Ds,imax + + Ds = 4 + Dv = 4 + imax = Ds/4 + + do i=1,imax + i1= 1+4*(i-1) + i2=i1+3 + + y1=sp(i1) + y2=sp(i1+1) + y3=sp(i1+2) + y4=sp(i1+3) + + x0(1,i)=y3 + x0(2,i)=y4 + x0(3,i)=y1 + x0(4,i)=y2 + + xx(1,i) = y4 + xx(2,i) = y3 + xx(3,i) = -y2 + xx(4,i) = -y1 + + xy(1,i)=(0d0,1d0)*y4 + xy(2,i)=-(0d0,1d0)*y3 + xy(3,i)=-(0d0,1d0)*y2 + xy(4,i)=(0d0,1d0)*y1 + + xz(1,i)=y3 + xz(2,i)=-y4 + xz(3,i)=-y1 + xz(4,i)=y2 + + x5(1,i)=y1 + x5(2,i)=y2 + x5(3,i)=-y3 + x5(4,i)=-y4 + enddo + + + do i=1,4 + spb2_Weyl(i)=v(1)*x0(i,1)-v(2)*xx(i,1)-v(3)*xy(i,1)-v(4)*xz(i,1) + enddo +end function + +function spi2_Weyl(v,sp) +implicit none +complex(8), intent(in) :: sp(:),v(:) +complex(8) :: spi2_Weyl(4) +complex(8) :: x0(4,4),xx(4,4),xy(4,4) +complex(8) :: xz(4,4),x5(4,4) +complex(8) :: y1,y2,y3,y4,bp,bm,cp,cm +integer :: i,i1,i2,i3,imax,Dv,Ds + + Ds = 4 + Dv = 4 + + imax = Ds/4 + + do i=1,imax + i1= 1+4*(i-1) + i2=i1+3 + + y1=sp(i1) + y2=sp(i1+1) + y3=sp(i1+2) + y4=sp(i1+3) + + x0(1,i)=y3 + x0(2,i)=y4 + x0(3,i)=y1 + x0(4,i)=y2 + + + xx(1,i) = -y4 + xx(2,i) = -y3 + xx(3,i) = y2 + xx(4,i) = y1 + + + xy(1,i)=(0d0,1d0)*y4 + xy(2,i)=-(0d0,1d0)*y3 + xy(3,i)=-(0d0,1d0)*y2 + xy(4,i)=(0d0,1d0)*y1 + + xz(1,i)=-y3 + xz(2,i)=y4 + xz(3,i)=y1 + xz(4,i)=-y2 + + x5(1,i)=y1 + x5(2,i)=y2 + x5(3,i)=-y3 + x5(4,i)=-y4 + enddo + do i=1,4 + spi2_Weyl(i)=v(1)*x0(i,1)-v(2)*xx(i,1) -v(3)*xy(i,1)-v(4)*xz(i,1) + enddo +end function + @@ -1034,187 +1191,187 @@ end subroutine convert_to_MCFM !might not be correct in this context !======================================================================== - SUBROUTINE EvaluateSpline(EvalPoint, SplineData, SplineDataLength, TheResult) - !SplineData: SplineDataLength by 2 array - ! SplineData(1:SplineDataLength,1) are the x values - ! SplineData(1:SplineDataLength,2) are the corresponding y values +SUBROUTINE EvaluateSpline(EvalPoint, SplineData, SplineDataLength, TheResult) +!SplineData: SplineDataLength by 2 array +! SplineData(1:SplineDataLength,1) are the x values +! SplineData(1:SplineDataLength,2) are the corresponding y values - IMPLICIT NONE +IMPLICIT NONE - INTEGER i,top,gdim,SplineDataLength - REAL(8) u,value,EvalPoint - REAL(8), intent(out) :: TheResult - REAL(8), dimension(SplineDataLength) :: bc,cc,dc - REAL(8) :: SplineData(1:SplineDataLength, 1:2) +INTEGER i,top,gdim,SplineDataLength +REAL(8) u,value,EvalPoint +REAL(8), intent(out) :: TheResult +REAL(8), dimension(SplineDataLength) :: bc,cc,dc +REAL(8) :: SplineData(1:SplineDataLength, 1:2) ! u value of M_H at which the spline is to be evaluated - gdim= SplineDataLength +gdim= SplineDataLength - CALL HTO_FMMsplineSingleHt(bc,cc,dc,top,gdim,SplineData(1:SplineDataLength,1),SplineData(1:SplineDataLength,2)) +CALL HTO_FMMsplineSingleHt(bc,cc,dc,top,gdim,SplineData(1:SplineDataLength,1),SplineData(1:SplineDataLength,2)) - u= EvalPoint - CALL HTO_Seval3SingleHt(u,bc,cc,dc,top,gdim,value,xc=SplineData(1:SplineDataLength,1),yc=SplineData(1:SplineDataLength,2)) +u= EvalPoint +CALL HTO_Seval3SingleHt(u,bc,cc,dc,top,gdim,value,xc=SplineData(1:SplineDataLength,1),yc=SplineData(1:SplineDataLength,2)) - TheResult= value +TheResult= value - RETURN +RETURN !----------------------------------------------------------------------- - CONTAINS +CONTAINS - SUBROUTINE HTO_FMMsplineSingleHt(b,c,d,top,gdim,xc,yc) +SUBROUTINE HTO_FMMsplineSingleHt(b,c,d,top,gdim,xc,yc) !--------------------------------------------------------------------------- - INTEGER k,n,i,top,gdim,l +INTEGER k,n,i,top,gdim,l - REAL(8), dimension(gdim) :: xc,yc - REAL(8), dimension(gdim) :: x,y +REAL(8), dimension(gdim) :: xc,yc +REAL(8), dimension(gdim) :: x,y - REAL(8), DIMENSION(gdim) :: b +REAL(8), DIMENSION(gdim) :: b ! linear coeff - REAL(8), DIMENSION(gdim) :: c +REAL(8), DIMENSION(gdim) :: c ! quadratic coeff. - REAL(8), DIMENSION(gdim) :: d +REAL(8), DIMENSION(gdim) :: d ! cubic coeff. - REAL(8) :: t - REAL(8),PARAMETER:: ZERO=0.0, TWO=2.0, THREE=3.0 +REAL(8) :: t +REAL(8),PARAMETER:: ZERO=0.0, TWO=2.0, THREE=3.0 ! The grid - n= gdim - FORALL(l=1:gdim) - x(l)= xc(l) - y(l)= yc(l) - ENDFORALL +n= gdim +FORALL(l=1:gdim) +x(l)= xc(l) +y(l)= yc(l) +ENDFORALL !.....Set up tridiagonal system......................................... ! b=diagonal, d=offdiagonal, c=right-hand side - d(1)= x(2)-x(1) - c(2)= (y(2)-y(1))/d(1) - DO k= 2,n-1 - d(k)= x(k+1)-x(k) - b(k)= TWO*(d(k-1)+d(k)) - c(k+1)= (y(k+1)-y(k))/d(k) - c(k)= c(k+1)-c(k) - END DO +d(1)= x(2)-x(1) +c(2)= (y(2)-y(1))/d(1) +DO k= 2,n-1 +d(k)= x(k+1)-x(k) +b(k)= TWO*(d(k-1)+d(k)) +c(k+1)= (y(k+1)-y(k))/d(k) +c(k)= c(k+1)-c(k) +END DO !.....End conditions. third derivatives at x(1) and x(n) obtained ! from divided differences....................................... - b(1)= -d(1) - b(n)= -d(n-1) - c(1)= ZERO - c(n)= ZERO - IF (n > 3) THEN - c(1)= c(3)/(x(4)-x(2))-c(2)/(x(3)-x(1)) - c(n)= c(n-1)/(x(n)-x(n-2))-c(n-2)/(x(n-1)-x(n-3)) - c(1)= c(1)*d(1)*d(1)/(x(4)-x(1)) - c(n)= -c(n)*d(n-1)*d(n-1)/(x(n)-x(n-3)) - END IF - - DO k=2,n ! forward elimination - t= d(k-1)/b(k-1) - b(k)= b(k)-t*d(k-1) - c(k)= c(k)-t*c(k-1) - END DO - - c(n)= c(n)/b(n) +b(1)= -d(1) +b(n)= -d(n-1) +c(1)= ZERO +c(n)= ZERO +IF (n > 3) THEN +c(1)= c(3)/(x(4)-x(2))-c(2)/(x(3)-x(1)) +c(n)= c(n-1)/(x(n)-x(n-2))-c(n-2)/(x(n-1)-x(n-3)) +c(1)= c(1)*d(1)*d(1)/(x(4)-x(1)) +c(n)= -c(n)*d(n-1)*d(n-1)/(x(n)-x(n-3)) +END IF + +DO k=2,n ! forward elimination +t= d(k-1)/b(k-1) +b(k)= b(k)-t*d(k-1) +c(k)= c(k)-t*c(k-1) +END DO + +c(n)= c(n)/b(n) ! back substitution ( makes c the sigma of text) - DO k=n-1,1,-1 - c(k)= (c(k)-d(k)*c(k+1))/b(k) - END DO +DO k=n-1,1,-1 +c(k)= (c(k)-d(k)*c(k+1))/b(k) +END DO !.....Compute polynomial coefficients................................... - b(n)= (y(n)-y(n-1))/d(n-1)+d(n-1)*(c(n-1)+c(n)+c(n)) - DO k=1,n-1 - b(k)= (y(k+1)-y(k))/d(k)-d(k)*(c(k+1)+c(k)+c(k)) - d(k)= (c(k+1)-c(k))/d(k) - c(k)= THREE*c(k) - END DO - c(n)= THREE*c(n) - d(n)= d(n-1) +b(n)= (y(n)-y(n-1))/d(n-1)+d(n-1)*(c(n-1)+c(n)+c(n)) +DO k=1,n-1 +b(k)= (y(k+1)-y(k))/d(k)-d(k)*(c(k+1)+c(k)+c(k)) +d(k)= (c(k+1)-c(k))/d(k) +c(k)= THREE*c(k) +END DO +c(n)= THREE*c(n) +d(n)= d(n-1) - RETURN +RETURN - END SUBROUTINE HTO_FMMsplineSingleHt +END SUBROUTINE HTO_FMMsplineSingleHt !------------------------------------------------------------------------ - SUBROUTINE HTO_Seval3SingleHt(u,b,c,d,top,gdim,f,fp,fpp,fppp,xc,yc) +SUBROUTINE HTO_Seval3SingleHt(u,b,c,d,top,gdim,f,fp,fpp,fppp,xc,yc) ! --------------------------------------------------------------------------- - REAL(8),INTENT(IN) :: u +REAL(8),INTENT(IN) :: u ! abscissa at which the spline is to be evaluated - INTEGER j,k,n,l,top,gdim +INTEGER j,k,n,l,top,gdim - REAL(8), dimension(gdim) :: xc,yc - REAL(8), dimension(gdim) :: x,y - REAL(8), DIMENSION(gdim) :: b,c,d +REAL(8), dimension(gdim) :: xc,yc +REAL(8), dimension(gdim) :: x,y +REAL(8), DIMENSION(gdim) :: b,c,d ! linear,quadratic,cubic coeff - REAL(8),INTENT(OUT),OPTIONAL:: f,fp,fpp,fppp +REAL(8),INTENT(OUT),OPTIONAL:: f,fp,fpp,fppp ! function, 1st,2nd,3rd deriv - INTEGER, SAVE :: i=1 - REAL(8) :: dx - REAL(8),PARAMETER:: TWO=2.0, THREE=3.0, SIX=6.0 +INTEGER, SAVE :: i=1 +REAL(8) :: dx +REAL(8),PARAMETER:: TWO=2.0, THREE=3.0, SIX=6.0 ! The grid - n= gdim - FORALL(l=1:gdim) - x(l)= xc(l) - y(l)= yc(l) - ENDFORALL +n= gdim +FORALL(l=1:gdim) +x(l)= xc(l) +y(l)= yc(l) +ENDFORALL !.....First check if u is in the same interval found on the ! last call to Seval............................................. - IF ( (i<1) .OR. (i >= n) ) i=1 - IF ( (u < x(i)) .OR. (u >= x(i+1)) ) THEN - i=1 +IF ( (i<1) .OR. (i >= n) ) i=1 +IF ( (u < x(i)) .OR. (u >= x(i+1)) ) THEN +i=1 ! binary search - j= n+1 - DO - k= (i+j)/2 - IF (u < x(k)) THEN - j= k - ELSE - i= k - ENDIF - IF (j <= i+1) EXIT - ENDDO - ENDIF +j= n+1 +DO +k= (i+j)/2 +IF (u < x(k)) THEN +j= k +ELSE +i= k +ENDIF +IF (j <= i+1) EXIT +ENDDO +ENDIF - dx= u-x(i) +dx= u-x(i) ! evaluate the spline - IF (Present(f)) f= y(i)+dx*(b(i)+dx*(c(i)+dx*d(i))) - IF (Present(fp)) fp= b(i)+dx*(TWO*c(i) + dx*THREE*d(i)) - IF (Present(fpp)) fpp= TWO*c(i) + dx*SIX*d(i) - IF (Present(fppp)) fppp= SIX*d(i) +IF (Present(f)) f= y(i)+dx*(b(i)+dx*(c(i)+dx*d(i))) +IF (Present(fp)) fp= b(i)+dx*(TWO*c(i) + dx*THREE*d(i)) +IF (Present(fpp)) fpp= TWO*c(i) + dx*SIX*d(i) +IF (Present(fppp)) fppp= SIX*d(i) - RETURN +RETURN - END SUBROUTINE HTO_Seval3SingleHt +END SUBROUTINE HTO_Seval3SingleHt - END SUBROUTINE EvaluateSpline +END SUBROUTINE EvaluateSpline function CenterWithStars(string, totallength) implicit none @@ -1234,5 +1391,37 @@ function CenterWithStars(string, totallength) end function +real function infinity() + implicit none + real :: x + x = 0 + infinity=-log(x) +end function infinity + +function CalculatesXsec(Process) + implicit none + integer :: Process + logical :: CalculatesXsec + if (Process.le.2) then + CalculatesXsec=.true. + elseif (Process.eq.50) then + CalculatesXsec=.false. + elseif (Process.eq.60) then + CalculatesXsec=.true. + elseif (Process.eq.61) then + CalculatesXsec=.true. + elseif (Process.eq.62) then + CalculatesXsec=.false. + elseif (Process.ge.66 .and. Process.le.69) then + CalculatesXsec=.true. + elseif (Process.eq.80 .or. Process.eq.90) then + CalculatesXsec=.false. + elseif (Process.ge.110 .and. Process.le.114) then + CalculatesXsec=.true. + else + print *, "Unknown process in CalculatesXsec", process + endif +end function CalculatesXsec + END MODULE diff --git a/MELA/fortran/mod_Parameters.F90 b/MELA/fortran/mod_Parameters.F90 index f3624e81..a64821f0 100755 --- a/MELA/fortran/mod_Parameters.F90 +++ b/MELA/fortran/mod_Parameters.F90 @@ -3,7 +3,7 @@ MODULE ModParameters save ! ! -character(len=6),parameter :: JHUGen_Version="v7.0.4" +character(len=6),parameter :: JHUGen_Version="v7.0.8" ! ! !===================================================== @@ -20,7 +20,7 @@ MODULE ModParameters integer, public :: FacScheme,RenScheme real(8), public :: MuFacMultiplier,MuRenMultiplier integer, public :: VegasIt1_default,VegasNc0_default,VegasNc1_default,VegasNc2_default -integer, public :: NumHistograms,RequestNLeptons,RequestOS,RequestOSSF +integer, public :: NumHistograms,RequestNLeptons,RequestOS,RequestOSSF,RequestNJets logical, public :: Unweighted,OffShellReson,OffShellV1,OffShellV2,ReadLHEFile,ConvertLHEFile,DoPrintPMZZ logical, public :: ReadCSmax,GenerateEvents,CountTauAsAny,HasLeptonFilter, FoundHiggsMass, FoundHiggsWidth integer, public :: WriteFailedEvents @@ -55,6 +55,12 @@ MODULE ModParameters integer, public :: LHAPDFMember, lenLHAPDFString ! lenLHAPDFString is needed in MELA integer, public :: PDFSet ! End PDFset variables +#if useCollier==1 +! COLLIER initialization variables +integer, public :: Collier_maxNLoopProps = -1 +integer, public :: Collier_maxRank = -1 +! End COLLIER initialization variables +#endif logical, public :: includeInterference, writegit real(8), public :: M_V,Ga_V, M_Vprime,Ga_Vprime real(8), public, parameter :: GeV=1d0/100d0 ! we are using units of 100GeV, i.e. Lambda=10 is 1TeV @@ -72,6 +78,7 @@ MODULE ModParameters integer, public :: Br_W_ud_counter=0 integer, public :: Br_counter(1:5,1:5)=0 integer, public :: LeptInEvent(0:8) = 0 +integer, public :: JetsInEvent(0:8) = 0 logical, public :: ReweightDecay = .false. integer, public :: UserSeed = 0 integer, public :: WidthScheme = 0 ! 1=running BW-width, 2=fixed BW-width (default), 3=Passarino's CPS @@ -127,12 +134,17 @@ MODULE ModParameters !===================================================== !cuts - should be set on the command line real(8), public :: pTjetcut = -1d0*GeV ! jet min pt, default is set in main (0 in VH, 15 GeV otherwise) +real(8), public :: etajetcut = -1d0 ! jet max |eta|, default is set in main (4 in offshell VBF, infinity elsewhere) +real(8), public :: detajetcut = -1d0 ! min difference in eta between jets (default 2 in VBF offshell, 0 elsewhere) real(8), public :: Rjet = -1d0 ! jet deltaR, anti-kt algorithm, default is set in main (0 in VH, 0.3 otherwise) real(8), public :: mJJcut = 0d0*GeV ! minimum mJJ for VBF, HJJ, bbH, VH real(8), public :: m4l_minmax(1:2) = (/ -1d0,-1d0 /)*GeV ! min and max for m_4l in off-shell VBF production; default is (-1,-1): m_4l ~ Higgs resonance (on-shell) logical, public :: includeGammaStar = .false. ! include offshell photons? logical, public :: includeVprime = .false. -real(8), public :: MPhotonCutoff = 4d0*GeV ! minimum |mass| for offshell photons when includeGammaStar = .true. +real(8), public :: MPhotonCutoff = -1d0*GeV ! minimum |mass_ll| for offshell photons when includeGammaStar = .true. or in VBF bkg +real(8), public :: pTlepcut = -1d0*GeV +real(8), public :: etalepcut = -1d0 +logical, public :: JetsOppositeEta = .true. !===================================================== @@ -1900,11 +1912,11 @@ FUNCTION IsAQuark(PartType) integer :: PartType IsAQuark=IsUpTypeQuark(PartType) .or. IsDownTypeQuark(PartType) END FUNCTION -FUNCTION IsLHEAQuark(PartType) +FUNCTION IsALHEQuark(PartType) implicit none -logical :: IsLHEAQuark +logical :: IsALHEQuark integer :: PartType - IsLHEAQuark=IsLHEUpTypeQuark(PartType) .or. IsLHEDownTypeQuark(PartType) + IsALHEQuark=IsLHEUpTypeQuark(PartType) .or. IsLHEDownTypeQuark(PartType) END FUNCTION FUNCTION IsALightQuark(PartType) @@ -1913,11 +1925,11 @@ FUNCTION IsALightQuark(PartType) integer :: PartType IsALightQuark=IsUpTypeLightQuark(PartType) .or. IsDownTypeQuark(PartType) END FUNCTION -FUNCTION IsLHEALightQuark(PartType) +FUNCTION IsALHELightQuark(PartType) implicit none -logical :: IsLHEALightQuark +logical :: IsALHELightQuark integer :: PartType - IsLHEALightQuark=IsLHEUpTypeLightQuark(PartType) .or. IsLHEDownTypeQuark(PartType) + IsALHELightQuark=IsLHEUpTypeLightQuark(PartType) .or. IsLHEDownTypeQuark(PartType) END FUNCTION @@ -1948,6 +1960,36 @@ FUNCTION IsALHELepton(PartType)! note that lepton means charged lepton here END FUNCTION +FUNCTION IsAGluon(PartType) +implicit none +logical :: IsAGluon +integer :: PartType + IsAGluon = (abs(PartType).eq.Glu_) +END FUNCTION + +FUNCTION IsALHEGluon(PartType) +implicit none +logical :: IsALHEGluon +integer :: PartType + IsALHEGluon = (abs(PartType).eq.convertLHE(Glu_)) +END FUNCTION + + +FUNCTION IsAJet(PartType) +implicit none +logical :: IsAJet +integer :: PartType + IsAJet = (IsALightQuark(PartType) .or. IsAGluon(PartType)) +END FUNCTION + +FUNCTION IsALHEJet(PartType) +implicit none +logical :: IsALHEJet +integer :: PartType + IsALHEJet = (IsALHELightQuark(PartType) .or. IsALHEGluon(PartType)) +END FUNCTION + + FUNCTION IsABoson(PartType) implicit none logical :: IsABoson @@ -2284,12 +2326,13 @@ subroutine ReadCommandLineArgument_logical(argument, argumentname, success, dest end subroutine ReadCommandLineArgument_logical -subroutine ReadCommandLineArgument_integer(argument, argumentname, success, dest, SetLastArgument, success2, success3, success4, success5) +subroutine ReadCommandLineArgument_integer(argument, argumentname, success, dest, SetLastArgument, success2, success3, success4, success5, multiply) implicit none character(len=*) :: argument, argumentname integer, intent(inout) :: dest logical, intent(inout) :: success logical, optional, intent(inout) :: SetLastArgument, success2, success3, success4, success5 +integer, optional, intent(in) :: multiply integer :: length if (present(SetLastArgument)) SetLastArgument=.false. @@ -2298,6 +2341,7 @@ subroutine ReadCommandLineArgument_integer(argument, argumentname, success, dest if( argument(1:length+1) .eq. trim(argumentname)//"=" ) then read(argument(length+2:len(argument)), *) dest + if (present(multiply)) dest = dest*multiply success=.true. if (present(SetLastArgument)) SetLastArgument=.true. if (present(success2)) success2=.true. @@ -2309,12 +2353,13 @@ subroutine ReadCommandLineArgument_integer(argument, argumentname, success, dest end subroutine ReadCommandLineArgument_integer -subroutine ReadCommandLineArgument_real8(argument, argumentname, success, dest, SetLastArgument, success2, success3, success4, success5) +subroutine ReadCommandLineArgument_real8(argument, argumentname, success, dest, SetLastArgument, success2, success3, success4, success5, multiply) implicit none character(len=*) :: argument, argumentname real(8), intent(inout) :: dest logical, intent(inout) :: success logical, optional, intent(inout) :: SetLastArgument, success2, success3, success4, success5 +real(8), optional, intent(in) :: multiply integer :: length if (present(SetLastArgument)) SetLastArgument=.false. @@ -2323,6 +2368,7 @@ subroutine ReadCommandLineArgument_real8(argument, argumentname, success, dest, if( argument(1:length+1) .eq. trim(argumentname)//"=" ) then read(argument(length+2:len(argument)), *) dest + if (present(multiply)) dest = dest*multiply success=.true. if (present(SetLastArgument)) SetLastArgument=.true. if (present(success2)) success2=.true. @@ -2334,13 +2380,15 @@ subroutine ReadCommandLineArgument_real8(argument, argumentname, success, dest, end subroutine ReadCommandLineArgument_real8 -subroutine ReadCommandLineArgument_complex8(argument, argumentname, success, dest, SetLastArgument, success2, success3, success4, success5) +subroutine ReadCommandLineArgument_complex8(argument, argumentname, success, dest, SetLastArgument, success2, success3, success4, success5, multiply, multiplyreal) implicit none character(len=*) :: argument, argumentname complex(8), intent(inout) :: dest real(8) :: re, im logical, intent(inout) :: success logical, optional, intent(inout) :: SetLastArgument, success2, success3, success4, success5 +complex(8), optional, intent(in) :: multiply +real(8), optional, intent(in) :: multiplyreal integer :: length if (present(SetLastArgument)) SetLastArgument=.false. @@ -2359,6 +2407,8 @@ subroutine ReadCommandLineArgument_complex8(argument, argumentname, success, des endif read(argument(length+2:len(argument)), *) re, im dest = dcmplx(re, im) + if (present(multiply)) dest = dest*multiply + if (present(multiplyreal)) dest = dest*multiplyreal success=.true. if (present(SetLastArgument)) SetLastArgument=.true. if (present(success2)) success2=.true. @@ -2795,6 +2845,28 @@ subroutine spinoru(p,za,zb,s) end subroutine spinoru !======================================================================== +!======================================================================== +! Common initialization functions that may be called multiple times if needed +! Check arXiv:1604.06792 for the parameters +subroutine InitCOLLIER(Nmax, Rmax) +#if useCollier==1 +use COLLIER +implicit none +integer, intent(in) :: Nmax, Rmax +integer :: supNmax, supRmax + supNmax = max(Nmax, Collier_maxNLoopProps) + supRmax = max(Rmax, Collier_maxRank) + if ((supNmax .gt. Collier_maxNLoopProps) .or. (supRmax .gt. Collier_maxRank)) then + call Init_cll(supNmax,supRmax,'') + call setMode_cll(1) + endif +#else +implicit none +integer, intent(in) :: Nmax, Rmax + return +#endif +end subroutine + END MODULE diff --git a/MELA/fortran/mod_VHiggs.F90 b/MELA/fortran/mod_VHiggs.F90 index f81ba7e0..dc48285c 100755 --- a/MELA/fortran/mod_VHiggs.F90 +++ b/MELA/fortran/mod_VHiggs.F90 @@ -7,7 +7,7 @@ module ModVHiggs !----- notation for subroutines - public :: EvalAmp_VHiggs,EvalUnpolAmpSq_gg_VH + public :: EvalAmp_VHiggs contains @@ -1710,221 +1710,5 @@ SUBROUTINE getSME(p,FermFlav,SME) - - - - ! p: g1, g2, Zl1, Zl2, H - SUBROUTINE EvalUnpolAmpSq_gg_VH(p,UnPolSqAmp) - use ModParameters - use ModMisc -#if useCollier==1 - use COLLIER -#endif - implicit none - real(8) :: p(1:4,1:4),UnPolSqAmp - complex(8) :: SP(1:7),SI(0:16),LC(1:15) - complex(8) :: qlI1,qlI2,qlI3,qlI4 - complex(8) :: ep1(1:4),ep2(1:4),cep3(1:4),p1(1:4),p2(1:4),p3(1:4) - complex(8) :: MasslessTri,MassiveTri,MassiveBox - real(8) :: shat,that,uhat,MT2,MH2,MZ2,MV2 - real(8) :: MuRen2,PreFac,PreFac2,IZ(3:4,9:9),smT,y,mTrans - integer :: eps,h1,h2,h3,i,j - complex(8) :: VVHg1,VVHg2,VVHg3 - complex(8) :: calc_MassiveBox,calc_MassiveBox_QP,calc_MassiveTensorBox - real(8),parameter :: E=dexp(1d0), ICol_sq=8d0 - -! real(8) :: sprod(1:4,1:4),p5hat(1:4),p6hat(1:4) -! complex(8) :: za(1:4,1:4),zb(1:4,1:4),a,b,c,d,e,f,v,w,x,y,z,o - - - -! prefactors and kinematics - PreFac = Pi**2 * (4d0*pi*alpha_QED) * (4d0*pi*alphas) /sitW/M_W *1d0/(16d0*pi**2) - PreFac2 = Pi**2 * dsqrt(4d0*pi*alpha_QED) * (4d0*pi*alphas) *m_Top**2/vev *1d0/(16d0*pi**2) - IZ(3,9) = 0.5d0*(aL_QUp-aR_QUp)/2d0/sitW/dsqrt(1d0-sitW**2) ! top - IZ(4,9) = -IZ(3,9) ! bottom - - p1(1:4) = p(1:4,1) - p2(1:4) = p(1:4,2) - p3(1:4) = p(1:4,3)+p(1:4,4) - MZ2 = M_Z**2 - MH2 = M_Reso**2 - MT2 = M_Top**2 - MV2 = get_MInv2(dble(p3(1:4))) - - VVHg1 = 1d0 - VVHg2 = 0d0/MH2 - VVHg3 = 0d0/MH2 - -! call ShiftMass2(p(1:4,3),p(1:4,4),M_Z,p5hat,p6hat)! projecting Z boson on-shell (because computation assumed p3^3=MZ^2) -! p3(1:4) = p5hat(1:4)+p6hat(1:4) - shat =+2d0*(p1.dot.p2) - that =-2d0*(p1.dot.p3)+MV2 - uhat =-2d0*(p2.dot.p3)+MV2 - mTrans = dsqrt( get_pT(dble(p3(1:4)))**2 + MV2 ) - smT = dsqrt(shat)*mTrans - y = Get_ETA(dble(p3(1:4))) - - -! print *, "check1",that,- smT* E**(-y) + MV2 -! print *, "check2",uhat,- smT* E**(+y) + MV2 -! pause - - -! p30(1:4) = p3(1:4) - (p3(1:4).dot.p3(1:4))/(2d0*(p3(1:4).dot.p1(1:4)))*p1(1:4) -! call my_spinoru(4,(/p(1:4,1),p(1:4,2),p5hat(1:4),p6hat(1:4)/),za,zb,sprod) -! a = za(1,2); b = za(1,3); c = za(1,4); d = za(2,3); e = za(2,4); f = za(3,4) -! v = zb(1,2); w = zb(1,3); x = zb(1,4); y = zb(2,3); z = zb(2,4); o = zb(3,4) - -! evaluate scalar integrals - SI(0) = 1d0 - eps = 0 -! ! check 1/eps^2 -! print *, "checking 1/eps^2" -! eps=-2 -! SI(0) = 0d0 - ! check 1/eps^1 -! print *, "checking 1/eps^1" -! eps=-1 -! SI(0) = 0d0 - - MuRen2 = Mu_Ren**2 -! print *, "checks" -! print *, PreFac,IZ,MZ2,MH2,MT2 -! print *, shat,that,uhat -! print *, p(1:4,1) -! print *, p(1:4,2) -! print *, p(1:4,3) -! print *, p(1:4,4) -! pause -#if useCollier==1 - call SetMuUV2_cll(MuRen2) - call SetMuIR2_cll(MuRen2) -#endif - - -! SI(1:16) = 0d0 -#if useQCDLoopLib==1 - SI(1) = qlI1(MT2,MuRen2,eps) - SI(2) = qlI2(zero, MT2,MT2,MuRen2,eps) - SI(3) = qlI2(MH2, MT2,MT2,MuRen2,eps) - SI(4) = qlI2(MV2, MT2,MT2,MuRen2,eps) - SI(5) = qlI2(shat, MT2,MT2,MuRen2,eps) - SI(6) = qlI2(that, MT2,MT2,MuRen2,eps) - SI(7) = qlI2(uhat, MT2,MT2,MuRen2,eps) - SI(8) = qlI3(zero,zero,shat, MT2,MT2,MT2,MuRen2,eps) - SI(9) = qlI3(zero,MV2,that, MT2,MT2,MT2,MuRen2,eps) - SI(10)= qlI3(zero,MV2,uhat, MT2,MT2,MT2,MuRen2,eps) - SI(11)= qlI3(zero,that,MH2, MT2,MT2,MT2,MuRen2,eps) - SI(12)= qlI3(zero,uhat,MH2, MT2,MT2,MT2,MuRen2,eps) - SI(13)= qlI3(MV2,shat,MH2, MT2,MT2,MT2,MuRen2,eps) - SI(14)= qlI4(zero,zero,MV2,MH2,shat,that, MT2,MT2,MT2,MT2,MuRen2,eps) - SI(15)= qlI4(zero,zero,MV2,MH2,shat,uhat, MT2,MT2,MT2,MT2,MuRen2,eps) - SI(16)= qlI4(zero,MV2,zero,MH2,that,uhat, MT2,MT2,MT2,MT2,MuRen2,eps) - -#else - - print *, "No loop integrals" -#endif - - -! helicity sum - UnPolSqAmp = 0d0 - do h1=-1,+1, 2 - do h2=-1,+1, 2 - do h3=-1,+1 !, 2 - -! polarization vectors - ep1(1:4) = pol_gluon_incoming(p(1:4,1),h1) - ep2(1:4) = pol_gluon_incoming(p(1:4,2),h2) -! cep3(1:4)= pol_Zff_outgoing(p(1:4,3),p(1:4,4),h3) * (-1d0)/( MV2-MZ2 + (0d0,1d0)*M_Z*Ga_Z ) -! if( h3.lt.0 ) then! CHECK THAT L-R ASSIGNMENTS ARE CORRECT !!! -! cep3(1:4) = cep3(1:4) * dsqrt(4d0*pi*alpha_QED)/2d0/sitW/dsqrt(1d0-sitW**2) * aR_lep -! else -! cep3(1:4) = cep3(1:4) * dsqrt(4d0*pi*alpha_QED)/2d0/sitW/dsqrt(1d0-sitW**2) * aL_lep -! endif - cep3(1:4) = pol_mass(p3(1:4),h3,outgoing=.true.) - - -! scalar products - SP(1) = (cep3.dot.ep1) - SP(2) = (cep3.dot.ep2) - SP(3) = (cep3.dot.p1) - SP(4) = (cep3.dot.p2) - SP(5) = ( ep1.dot.ep2) - SP(6) = ( ep1.dot.p3) - SP(7) = ( ep2.dot.p3) - -! Levi-Civita tensors - LC(1) = LeviCiv(p1,ep1,ep2,cep3) - LC(2) = LeviCiv(p1,p2,ep1,cep3) - LC(3) = LeviCiv(p1,p2,ep1,ep2) - LC(4) = LeviCiv(p1,p2,ep2,cep3) - LC(5) = LeviCiv(p1,p2,p3,cep3) - LC(6) = LeviCiv(p1,p2,p3,ep1) - LC(7) = LeviCiv(p1,p2,p3,ep2) - LC(8) = LeviCiv(p1,p3,ep1,cep3) - LC(9) = LeviCiv(p1,p3,ep1,ep2) - LC(10)= LeviCiv(p1,p3,ep2,cep3) - LC(11)= LeviCiv(p2,ep1,ep2,cep3) - LC(12)= LeviCiv(p2,p3,ep1,cep3) - LC(13)= LeviCiv(p2,p3,ep1,ep2) - LC(14)= LeviCiv(p2,p3,ep2,cep3) - LC(15)= LeviCiv(p3,ep1,ep2,cep3) - - - - -! Loop amplitudes - MasslessTri = (2*MZ2*SI(0)*(-4*VVHg1*LC(3)*(SP(3) + SP(4)) + 4*smT*VVHg2*Cosh(y)*LC(3)*(SP(3) + SP(4)) & - + shat*(-(VVHg1*(LC(1) + 5*LC(11))) + VVHg2*(LC(9) + 5*LC(13))*(SP(3) + SP(4))) & - + 3*shat**2*VVHg3*(SP(2)*SP(6) - SP(1)*SP(7))))/(3*shat*(-MH2 + MV2 + MZ2 - 2*smT*Cosh(y))) -! - MassiveTri = (-2*(2*smT*VVHg2*LC(3)*(MZ2*(shat*SI(0) - SI(1) + MT2*(SI(0) - 2*SI(2) + 3*SI(5))) + & - 3*MT2*(MZ2 - shat)*shat*SI(8))*(SP(3) + SP(4)) + 2*E**(2*y)*smT*VVHg2*LC(3)* & - (MZ2*(shat*SI(0) - SI(1) + MT2*(SI(0) - 2*SI(2) + 3*SI(5))) + 3*MT2*(MZ2 - shat)*shat*SI(8))*(SP(3) + SP(4)) + & - E**y*(-2*MT2*(-6*shat**2*(VVHg1 + shat*VVHg2)*LC(3)*SI(8)*(SP(3) + SP(4)) + MZ2*(2*VVHg1*LC(3)*(SI(0) - 2*SI(2) + 3*SI(5))* & - (SP(3) + SP(4)) + 6*shat**2*VVHg2*LC(3)*SI(8)*(SP(3) + SP(4)) + shat*(VVHg2*(LC(9) - LC(13))*(SI(0) - 2*SI(2) + 3*SI(5))* & - (SP(3) + SP(4)) + VVHg1*(-((LC(1) - LC(11))*(SI(0) - 2*SI(2) + 3*SI(5))) + & - 6*LC(3)*SI(8)*(SP(3) + SP(4)))))) + MZ2*(4*VVHg1*LC(3)*SI(1)*(SP(3) + SP(4)) + & - shat**2*SI(0)*(-(VVHg1*(LC(1) + 5*LC(11))) + VVHg2*(LC(9) + 5*LC(13))*(SP(3) + SP(4))) - & - 2*shat*(-(VVHg2*(LC(9) - LC(13))*SI(1)*(SP(3) + SP(4))) + VVHg1*((LC(1) - LC(11))*SI(1) + & - 2*LC(3)*SI(0)*(SP(3) + SP(4)))) + 3*shat**3*VVHg3*SI(0)*(SP(2)*SP(6) - SP(1)*SP(7))))))/(3*E**y*shat**2*(MH2 - MV2 - MZ2 + 2*smT*Cosh(y))) - -! MassiveBox = calc_MassiveBox(MV2,MH2,MT2,shat,smT,y,LC,SP,SI,kappa,kappa_tilde) -! MassiveBox = calc_MassiveBox_QP(MV2,MH2,MT2,shat,smT,y,LC,SP,SI,kappa,kappa_tilde) - MassiveBox = calc_MassiveTensorBox(MV2,MT2,MH2,shat,that,uhat,smT,y,LC,SP,SI,kappa,kappa_tilde) - - - MasslessTri = MasslessTri * PreFac * IZ(4,9) - MassiveTri = MassiveTri * PreFac * IZ(3,9) - MassiveBox = MassiveBox * PreFac2* IZ(3,9) - - UnPolSqAmp = UnPolSqAmp + ICol_sq * cdabs( MasslessTri*1 + MassiveTri*1 + MassiveBox*1 )**2 - - - -! print *, h1,h2,h3 -! print *, "MasslessTri",MasslessTri -! print *, "MassiveTri",MassiveTri -! print *, "MassiveBox",MassiveBox -! pause - - enddo - enddo - enddo - - - RETURN - END SUBROUTINE - - - - - - - - - end module ModVHiggs !!--YaofuZhou----------------------------------------- diff --git a/MELA/setup.sh b/MELA/setup.sh index 736bae63..8d9fc208 100755 --- a/MELA/setup.sh +++ b/MELA/setup.sh @@ -4,6 +4,9 @@ set -euo pipefail cd $(dirname ${BASH_SOURCE[0]}) eval $(scram ru -sh) + +bash COLLIER/setup.sh "$@" + if [[ $# > 0 ]] && [[ "$1" == *"clean"* ]];then scramv1 b "$@" pushd ${CMSSW_BASE}/src/ZZMatrixElement/MELA/fortran/ diff --git a/MELA/src/TNumericUtil.cc b/MELA/src/TNumericUtil.cc index c1cf7bdc..b9c4df48 100755 --- a/MELA/src/TNumericUtil.cc +++ b/MELA/src/TNumericUtil.cc @@ -4,32 +4,32 @@ using namespace std; -void TNumericUtil::PermutationGenerator(int n, int k, std::vector>& perms, int valbegin, int increment){ - if (n<=0 || k<=0 || k>n) return; - vector d(n); - iota(d.begin(), d.end(), 1); - do{ - vector aperm; - for (int i=0; i>& perms, int valbegin, int increment){ + if (n<=0 || k<=0 || k>n) return; + vector d(n); + iota(d.begin(), d.end(), 1); + do{ + vector aperm; + for (int i=0; i>& perms, int valbegin, int increment){ - if (n==k && k>0){ - vector d(n); - iota(d.begin(), d.end(), 1); - perms.push_back(d); - } - else{ - PermutationGenerator(n, k, perms); - for (auto& perm:perms) sort(perm.begin(), perm.end()); - for (auto ip=perms.begin(); ip!=perms.end(); ip++){ - for (auto jp=perms.rbegin(); (jp.base()-1)!=ip; jp++){ - if (*(ip)==*(jp)) perms.erase(jp.base()-1); - } - } - } - if (valbegin!=1 || increment!=1){ for (auto& p:perms){ for (auto& pp:p) pp = increment*(pp-1)+valbegin; } } +void TNumericUtil::CombinationGenerator(int n, int k, std::vector>& perms, int valbegin, int increment){ + if (n==k && k>0){ + vector d(n); + iota(d.begin(), d.end(), 1); + perms.push_back(d); + } + else{ + PermutationGenerator(n, k, perms); + for (auto& perm:perms) sort(perm.begin(), perm.end()); + for (auto ip=perms.begin(); ip!=perms.end(); ip++){ + for (auto jp=perms.rbegin(); (jp.base()-1)!=ip; jp++){ + if (*(ip)==*(jp)) perms.erase(jp.base()-1); + } + } + } + if (valbegin!=1 || increment!=1){ for (auto& p:perms){ for (auto& pp:p) pp = increment*(pp-1)+valbegin; } } } diff --git a/MELA/test/loadMELA.C b/MELA/test/loadMELA.C index 21419f29..02d63977 100755 --- a/MELA/test/loadMELA.C +++ b/MELA/test/loadMELA.C @@ -1,9 +1,12 @@ { + TString LIBCOLLIER="libcollier.so"; TString LIBMCFM="libmcfm_704.so"; + gSystem->AddIncludePath("-I$ROOFITSYS/include/"); gSystem->AddIncludePath("-I$CMSSW_BASE/interface/"); gSystem->AddIncludePath("-I$CMSSW_BASE/src/"); gSystem->AddIncludePath("-I$CMSSW_BASE/src/ZZMatrixElement/MELA/interface/"); + gSystem->Load("$CMSSW_BASE/src/ZZMatrixElement/MELA/data/$SCRAM_ARCH/" + LIBCOLLIER); gSystem->Load("libZZMatrixElementMELA.so"); gSystem->Load("$CMSSW_BASE/src/ZZMatrixElement/MELA/data/$SCRAM_ARCH/" + LIBMCFM); }