subroutine upinit C...Double precision and integer declarations. IMPLICIT DOUBLE PRECISION(A-H, O-Z) IMPLICIT INTEGER(I-N) C...Commonblocks. C...The event record. COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5) C...Parameters. COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) C...Particle properties + some flavour parameters. COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4) C...User process initialization commonblock. INTEGER MAXPUP PARAMETER (MAXPUP=100) INTEGER IDBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,LPRUP DOUBLE PRECISION EBMUP,XSECUP,XERRUP,XMAXUP COMMON/HEPRUP/IDBMUP(2),EBMUP(2),PDFGUP(2),PDFSUP(2), &IDWTUP,NPRUP,XSECUP(MAXPUP),XERRUP(MAXPUP),XMAXUP(MAXPUP), &LPRUP(MAXPUP) C...Set default incoming beams, protons, and energy, 7 TeV. data idbmup/2212,2212/ data ebmup/0.7d4,0.7d4/ SAVE /HEPRUP/ double precision tau_c,y_c,z_c common/gen_c/tau_c(3),y_c(3),z_c(3) character*80 higgsfile common/hf/higgsfile C..Setup the number of processes and the type nprup=2 lprup(1)=101 lprup(2)=102 idwtup=1 C...A low value of the maximul differential cross-section is used as defult xmaxup(1)=1e-1 xmaxup(2)=1e-1 C...Some coefficients for phases space generation tau_c(1)=0.25d0 tau_c(2)=0.0d0 tau_c(3)=0.75d0 y_c(1)=0d0 y_c(2)=1d0 y_c(3)=0d0 z_c(1)=0d0 z_c(2)=0d0 z_c(3)=1d0 return end subroutine upcoup() C...Double precision and integer declarations. IMPLICIT DOUBLE PRECISION(A-H, O-Z) IMPLICIT INTEGER(I-N) C...Parameters. COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) C...Particle properties + some flavour parameters. COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4) C...Supersymmetry parameters. COMMON/PYMSSM/IMSS(0:99),RMSS(0:99) C...User process initialization commonblock. INTEGER MAXPUP PARAMETER (MAXPUP=100) INTEGER IDBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,LPRUP DOUBLE PRECISION EBMUP,XSECUP,XERRUP,XMAXUP COMMON/HEPRUP/IDBMUP(2),EBMUP(2),PDFGUP(2),PDFSUP(2), &IDWTUP,NPRUP,XSECUP(MAXPUP),XERRUP(MAXPUP),XMAXUP(MAXPUP), &LPRUP(MAXPUP) SAVE /HEPRUP/ integer ipybbwh double precision o_m common/pybbwh/o_m(3,3),ipybbwh data ipybbwh/0/ complex*16 ghichw,ghibb common/pybbwhg/ghichw(3),ghibb(3) double precision alpha,beta double precision s,tau,y,z,dsigma,ps double precision dmHc C...Extrat alpha and beta from Pythia common block ALPHA=RMSS(18) BETA=atan(RMSS(5)) if (ipybbwh.eq.0) then C...Set up the Higgs mixing matrix o_m(1,1)=-dsin(alpha) o_m(1,2)= dcos(alpha) o_m(1,3)=0d0 o_m(2,1)=o_m(2,1) o_m(2,2)=dsin(alpha) o_m(2,3)=0d0 o_m(3,1)=0d0 o_m(3,2)=0d0 o_m(3,3)=1 C...Some consistency checks if (dabs((-o_m(1,1)/dcos(beta))/paru(161)-1).gt.0.1) then write(*,*) $ 'Warning: RMSS(18) is inconsistent with h0 coupling' write(*,*) $ (-o_m(1,1)/dcos(beta)),paru(161) endif if (dabs((-o_m(1,2)/dcos(beta))/paru(171)-1).gt.0.1) then write(*,*) $ 'Warning: RMSS(18) is inconsistent with H0 coupling' write(*,*) $ (-o_m(1,2)/dcos(beta)),paru(171) endif if (dabs((dtan(beta))/paru(181)-1).gt.0.1) then write(*,*) $ 'Warning: RMSS(18) is inconsistent with A0 coupling' write(*,*) (dtan(beta)),paru(181) endif C...Set some Higgs couplings not normaly set paru(197)=dsin(alpha-beta) paru(198)=1 endif C...Calculate the couplings used in the matrix element do 10 i=1,3 ghichw(i)=dcmplx(o_m(2,i)*dcos(beta)-o_m(1,i)*dsin(beta), $ o_m(3,i)) ghibb(i)=dcmplx(o_m(1,i),o_m(3,i)*dsin(beta)) 10 continue C...Calculate xmaxup from off-shell resonant part s=(ebmup(1)+ebmup(2))**2 ps=1d0 dmHc=(pmas(37,1)-pmas(24,1))/2d0 if (dmHc.lt.pmas(5,1)) dmHc=pmas(37,1) call tau_gen(0.5d0,0.5d0,(PMAS(24,1)+dmHc)**2/s,1d0,tau,psfac) ps=ps*psfac call y_gen(0.5d0,0.5d0,log(tau)/2d0,-log(tau)/2d0,y,psfac) ps=ps*psfac call z_gen(0.5d0,0.5d0,-1d0,1d0,PMAS(24,1),dmHc,s,tau,z,psfac) ps=ps*psfac xwgtup=ps*dsigma(1,s,idbmup(1),idbmup(2), $ tau,y,z,PMAS(24,1),dmHc)*2d12 xmaxup(1)=xwgtup xmaxup(2)=xwgtup return end subroutine upevnt C...Double precision and integer declarations. IMPLICIT DOUBLE PRECISION(A-H, O-Z) IMPLICIT INTEGER(I-N) C...Commonblocks. C...The event record. COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5) C...Parameters. COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) C...Particle properties + some flavour parameters. COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4) C...User process event common block. INTEGER MAXNUP PARAMETER (MAXNUP=500) INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,IDUP(MAXNUP), &ISTUP(MAXNUP),MOTHUP(2,MAXNUP),ICOLUP(2,MAXNUP),PUP(5,MAXNUP), &VTIMUP(MAXNUP),SPINUP(MAXNUP) SAVE /HEPEUP/ C...User process initialization commonblock. INTEGER MAXPUP PARAMETER (MAXPUP=100) INTEGER IDBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,LPRUP DOUBLE PRECISION EBMUP,XSECUP,XERRUP,XMAXUP COMMON/HEPRUP/IDBMUP(2),EBMUP(2),PDFGUP(2),PDFSUP(2), &IDWTUP,NPRUP,XSECUP(MAXPUP),XERRUP(MAXPUP),XMAXUP(MAXPUP), &LPRUP(MAXPUP) double precision s,tau,y,z,dsigma,ps double precision tau_min,tau_max,y_min,y_max,z_min,z_max double precision shat,beta34,pT,BW_dmW,BW_dmHc C...Mass-dependet widths double precision dmW_min,dmW_max,dmHc_min,dmHc_max double precision gHpmin,gHpmax,gammaHp common/gammafh/gHpmin,gHpmax,gammaHp(100) logical first data first/.true./ save first C...If this is the first call setup up the couplings if (first) then first=.false. call upcoup() endif C...Set up outgoing event if (idwtup.eq.4) then if (pyr(0).lt.0.5d0) then idprup=101 else idprup=102 endif endif nup=4 do 10 i=1,2 istup(i)=-1 mothup(1,i)=0 mothup(2,i)=0 10 continue do 20 i=3,4 istup(i)=1 mothup(1,i)=1 mothup(2,i)=2 icolup(1,i)=0 icolup(2,i)=0 20 continue do 30 i=1,4 spinup(i)=9 vtimup(i)=0 30 continue aqedup=-1 aqcdup=-1 C...Random choice of where the b-quark comes from iqqbar=0 if (pyr(0).lt.0.5d0) iqqbar=1 if (iqqbar.eq.0) then idup(1)=5 idup(2)=-5 icolup(1,1)=1 icolup(1,2)=0 icolup(2,1)=0 icolup(2,2)=1 else idup(1)=-5 idup(2)=5 icolup(1,1)=0 icolup(1,2)=1 icolup(2,1)=1 icolup(2,2)=0 endif C...Set outgoing charges from selected subprocess if (idprup-100.eq.1) then idup(3)=-24 idup(4)=37 else idup(3)=24 idup(4)=-37 endif C...Start of phases space generation s=(ebmup(1)+ebmup(2))**2 ps=1d0 C...Generate W and H+ masses call mass_W_Hc(s,pyr(0),pyr(0),pyr(0),pyr(0),BW_dmW,BW_dmHc,psfac) ps=ps*psfac if (BW_dmW.lt.0d0) then xwgtup=0d0 return endif C...Generate tau tau_min=(BW_dmW+BW_dmHc)**2/s tau_max=1d0 if (tau_max.lt.tau_min) then xwgtup=0d0 return endif call tau_gen(pyr(0),pyr(0),tau_min,tau_max,tau,psfac) ps=ps*psfac C...Generate y y_min=log(tau)/2d0 y_max=-log(tau)/2d0 call y_gen(pyr(0),pyr(0),y_min,y_max,y,psfac) ps=ps*psfac x1=exp(y)*sqrt(tau) x2=sqrt(tau)/exp(y) if (x1.lt.0d0.or.x2.gt.1d0.or.x2.lt.0d0.or.x2.gt.1d0) then write(*,*) 'Something went wrong' write(*,*) 'dmW, mH, x1, x2:',BW_dmW,BW_dmHc,x1,x2 write(*,*) 'dmW_min, dmW_max, mH_min, mH_max, sqrt(s):' $ ,dmW_min,dmW_max,dmHc_min,dmHc_max,sqrt(s) write(*,*) 'tau, tau_min, tau_max: ',tau,tau_min,tau_max write(*,*) 'y, y_min, y_max: ',y,y_min,y_max endif C...Generate z z_min=-1d0 z_max=1d0 call z_gen(pyr(0),pyr(0),z_min,z_max,BW_dmW,BW_dmHc,s,tau,z,psfac) ps=ps*psfac C...Calculate the cross-section at the phase space point xwgtup=ps*dsigma(idprup-100+2*iqqbar,s,idbmup(1),idbmup(2), $ tau,y,z,BW_dmW,BW_dmHc)*2d12 C...Make the b-quark direction correct if (iqqbar.eq.1) z=-z C...Set up the outgoing partons shat=s*tau beta34=sqrt((1d0-BW_dmW**2/shat-BW_dmHc**2/shat)**2 $ -4d0*BW_dmW**2*BW_dmHc**2/shat**2) that=-1d0/2d0*((shat-BW_dmW**2-BW_dmHc**2)-shat*beta34*z) pT=sqrt(shat)/2d0*beta34*sqrt(1d0-z**2) Ebeam=sqrt(s/4d0) k(101,1)=1 p(101,1)=0d0 p(101,2)=0d0 p(101,3)=sqrt(shat)/2d0 p(101,4)=sqrt(shat)/2d0 p(101,5)=0d0 k(102,1)=1 p(102,1)=0d0 p(102,2)=0d0 p(102,3)=-sqrt(shat)/2d0 p(102,4)=sqrt(shat)/2d0 p(102,5)=0d0 k(103,1)=1 p(103,1)=pT p(103,2)=0d0 p(103,3)=sqrt(shat)/2d0*beta34*z p(103,4)=(shat+(BW_dmW**2-BW_dmHc**2))/2d0/sqrt(shat) p(103,5)=BW_dmW k(104,1)=1 p(104,1)=-pT p(104,2)=0d0 p(104,3)=-sqrt(shat)/2d0*beta34*z p(104,4)=(shat-(BW_dmW**2-BW_dmHc**2))/2d0/sqrt(shat) p(104,5)=BW_dmHc C...Boost to the correct frame n=104 call pyrobo(101,104,0d0,0d0,0d0,0d0, $ (exp(2d0*y)-1d0)/(exp(2d0*y)+1)) call pyrobo(101,104,0d0,pyr(0)*paru(2),0d0,0d0,0d0) C...Fill the outgoing event do 100 i=1,4 do 100 j=1,5 100 pup(j,i)=p(100+i,j) return end subroutine mass_W_Hc(s,rand1,rand2,rand3,rand4,BW_dmW,BW_dmHc,ps) C...Double precision and integer declarations. IMPLICIT DOUBLE PRECISION(A-H, O-Z) IMPLICIT INTEGER(I-N) C...Commonblocks. C...Particle properties + some flavour parameters. COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4) double precision dmW,dmHc C...Mass-dependet widths dimension wdtp(0:400),wdte(0:400,0:5) double precision dmW_min,dmW_max,dmHc_min,dmHc_max double precision gHpmin,gHpmax,gammaHp common/gammafh/gHpmin,gHpmax,gammaHp(100) C...Fixed masses dmW=pmas(24,1) dmHc=pmas(37,1) C...BW mass for W dmW=pmas(24,1) GW=pmas(24,2) dmW_min=0d0 dmW_max=sqrt(s) a_low=atan((dmW_min**2-dmW**2)/(dmW*GW)) a_high=atan((dmW_max**2-dmW**2)/(dmW*GW)) c_W=a_high-a_low BW_dmW=sqrt(dmW**2+dmW*GW*tan(a_low+rand1*(a_high-a_low))) C...BW mass for Hc dmHc=pmas(37,1) GHc=pmas(37,2) dmHc_min=0d0 dmHc_max=sqrt(s)-BW_dmW a_low=atan((dmHc_min**2-dmHc**2)/(dmHc*GHc)) a_high=atan((dmHc_max**2-dmHc**2)/(dmHc*GHc)) c_Hc=a_high-a_low BW_dmHc=sqrt(dmHc**2+dmHc*GHc*tan(a_low+rand2*(a_high-a_low))) C...Extra term for off-shell resonant contribution Res_min=pmas(6,1)*2d0/6d0 if (Res_min.lt.dmHc_min) Res_min=dmHc_min Res_max=(dmHc-dmW) if (Res_max.lt.Res_min) then c_R_Hc=0d0 Res_max=Res_min-1d0 else BW_mR=Res_min+rand3*(Res_max-Res_min) c_R_Hc=1d0/(pmas(37,1)*pmas(37,2))/(a_high-a_low) c_R_Hc=c_R_Hc/(c_R_Hc+10d0/2d0/(Res_max-Res_min)**2) if (rand4.lt.c_R_Hc) BW_dmHc=BW_mR endif C...Phase space factor ps=1d0 C...Mass-dependent width for W given by pywidt for phase spacec factor GMW=pmas(24,1)*pmas(24,2) hBWW=GMW/((BW_dmW**2-dmW**2)**2+GMW**2) call pywidt(24,BW_dmW**2,wdtp,wdte) GMW_c=BW_dmW*wdtp(0) hBWW_c=GMW_c/((BW_dmW**2-dmW**2)**2+GMW_c**2) ps=ps*hBWW_c/hBWW C..Mass-dependent width for Hc given by Hpwid GMHc=pmas(37,1)*pmas(37,2) hBWHc=GMHc/((BW_dmHc**2-dmHc**2)**2+GMHc**2) GMHc_c=BW_dmHc*Hpwid(BW_dmHc) hBWHc_c=GMHc_c/((BW_dmHc**2-dmHc**2)**2+GMHc_c**2) C...Phase space factor for mass coming from BW or off-shell resonant part if (BW_dmHc.gt.Res_min.and.BW_dmHc.lt.Res_max) then ps=ps*hBWHc_c/(a_high-a_low)/ $ ((1d0-c_R_Hc)*hBWHc/(a_high-a_low) $ +c_R_Hc/2d0/BW_dmHc/(Res_max-Res_min)) else ps=ps*hBWHc_c/((1d0-c_R_Hc)*hBWHc) endif return end double precision function Hpwid(MHp) C...Double precision and integer declarations. IMPLICIT DOUBLE PRECISION(A-H, O-Z) IMPLICIT INTEGER(I-N) C...Particle properties + some flavour parameters. COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4) double precision MHp,wdtp,wdte dimension wdtp(0:400),wdte(0:400,0:5) integer ipybbwh double precision o_m common/pybbwh/o_m(3,3),ipybbwh if (ipybbwh.eq.0) then call pywidt(37,MHp**2,wdtp,wdte) Hpwid=wdtp(0) else Hpwid=PMAS(37,2) endif return end subroutine tau_gen(rand1,rand2,tau_min,tau_max,tau,psfac) double precision rand1,rand2,tau_min,tau_max,tau,psfac double precision a double precision tau_c,y_c,z_c common/gen_c/tau_c(3),y_c(3),z_c(3) a=-3d0 C...tau from 1/tau if (rand1.lt.tau_c(1)) then tau=exp(log(tau_min)+rand2*(log(tau_max/tau_min))) C...tau from 1/tau**2 elseif (rand1.lt.(tau_c(1)+tau_c(2))) then tau=1d0/(1d0/tau_min+rand2*(1d0/tau_max-1d0/tau_min)) C...tau from tau**a else tau=(tau_min**(a+1d0)+ $ rand2*(tau_max**(a+1d0)-tau_min**(a+1d0)))**(1d0/(a+1d0)) endif C...phase space factor psfac=1d0/( $ 1d0/tau/log(tau_max/tau_min)*tau_c(1)+ $ 1d0/tau**2/(1d0/tau_min-1d0/tau_max)*tau_c(2)+ $ tau**a/(tau_max**(a+1d0)-tau_min**(a+1d0))*(a+1d0)*tau_c(3)) return end subroutine y_gen(rand1,rand2,y_min,y_max,y,psfac) double precision rand1,rand2,y_min,y_max,y,psfac double precision tau_c,y_c,z_c common/gen_c/tau_c(3),y_c(3),z_c(3) C...y from flat if (rand1.lt.y_c(1)) then y=y_min+(y_max-y_min)*rand2 C...y from else y=log(tan(atan(exp(y_min))+ $ rand2*(atan(exp(y_max))-atan(exp(y_min))))) endif C...phase spcae factor psfac=1d0/( $ 1d0/(y_max-y_min)*y_c(1)+ $ 1d0/(cosh(y)*2d0*(atan(exp(y_max))-atan(exp(y_min))))*y_c(2)) return end subroutine z_gen(rand1,rand2,z_min,z_max,dmW,dmHc,s,tau,z,psfac) double precision rand1,rand2,z_min,z_max,dmW,dmHc,s,tau,z,psfac double precision a,ce double precision tau_c,y_c,z_c common/gen_c/tau_c(3),y_c(3),z_c(3) ce=2d0 a=1d0+2d0*dmW**2*dmHc**2/s**2/tau**2 C...z from flat if (rand1.lt.z_c(1)) then z=z_min+(z_max-z_min)*rand2 C...z from elseif (rand1.lt.(z_c(1)+z_c(2))) then z=a-exp(log(a-z_min)+rand2*(log(a-z_max)-log(a-z_min))) C...z from else z=log(exp(ce*z_min)+rand2*(exp(ce*z_max)-exp(ce*z_min)))/ce endif C...phase spcae factor psfac=1d0/(1d0/(z_max-z_min)*z_c(1)+ $ 1d0/(a-z)/(-log(a-z_max)+log(a-z_min))*z_c(2)+ $ (exp(ce*z))/(exp(ce*z_max)-exp(ce*z_min))*ce*z_c(3)) return end function dsigma(proc,s,p1,p2,tau,y,z,dmW,dmHc) C...Double precision and integer declarations. IMPLICIT DOUBLE PRECISION(A-H, O-Z) IMPLICIT INTEGER(I-N) C...Parameters. COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) integer proc,p1,p2 double precision s,tau,y,z,dsigma double precision x1,x2,shat,that,beta34,pT double precision dmW,dmHc,me_bbWmHp,me_bbWpHm double precision xpq1,xpq2,Q2 C...PDF results dimension xpq1(-25:25),xpq2(-25:25) dsigma=1d0 C...phase space variables x1=exp(y)*sqrt(tau) x2=sqrt(tau)/exp(y) shat=s*x1*x2 beta34=sqrt((1d0-dmW**2/shat-dmHc**2/shat)**2 $ -4d0*dmW**2*dmHc**2/shat**2) that=-1d0/2d0*((shat-dmW**2-dmHc**2)-shat*beta34*z) uhat=-1d0/2d0*((shat-dmW**2-dmHc**2)+shat*beta34*z) pT=sqrt(shat)/2d0*beta34*sqrt(1d0-z**2) C...PDF Q2=((dmW+dmHc)/2.)**2 SCALUP=sqrt(Q2) call pypdfu(p1,x1,Q2,xpq1) call pypdfu(p2,x2,Q2,xpq2) C...phase space factor if (proc.gt.2) then dsigma=dsigma/tau*shat*beta34/2d0 $ *xpq1(-5)*xpq2(5) else dsigma=dsigma/tau*shat*beta34/2d0 $ *xpq1(5)*xpq2(-5) endif C...Matrix element for the two processes if (mod(proc,2).eq.1) then dsigma=dsigma*me_bbWmHp(shat,that,pT,dmW,dmHc) else dsigma=dsigma*me_bbWpHm(shat,that,pT,dmW,dmHc) endif dsigma=dsigma*1d-3*paru(5) return end double precision function me_bbWmHp(s,t,pT,dmW,mHp) C...Double precision and integer declarations. IMPLICIT DOUBLE PRECISION(A-H, O-Z) IMPLICIT INTEGER(I-N) C...Parameters. COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) C...Particle properties + some flavour parameters. COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4) C...Supersymmetry parameters. COMMON/PYMSSM/IMSS(0:99),RMSS(0:99) double precision s,t,pT,lambda,beta,ss,ts double precision mbrun,mHp,dmW,mt,mtrun,Gf,pi,Q2 complex*16 ghichw,ghibb common/pybbwhg/ghichw(3),ghibb(3) complex*16 prop integer hkf(1:3) data hkf/25,35,36/ C...Propagator prop(s,kf)=(1d0/dcmplx( $ s-pmas(kf,1)**2, $ pmas(kf,1)*pmas(kf,2))) lambda(x,y,z)=x**2+y**2+z**2-2*(x*y+y*z+z*x) beta=atan(RMSS(5)) Q2=(dmW+mHp)**2 mbrun=pymrun(5,Q2) mt=pmas(6,1) mtrun=pymrun(6,Q2) Gf=paru(1)*pyalem(Q2)/(sqrt(2d0)*paru(102)*pmas(24,1)**2) pi=paru(1) C...s-channel term ss=0d0 dss=0d0 do 20 i=1,3 do 10 j=1,3 ss=ss+dreal( $ ghichw(i)*dconjg(ghichw(j))* $ prop(s,hkf(i))*dconjg(prop(s,hkf(j)))* $ dreal(ghibb(i)*dconjg(ghibb(j))))/dcos(beta)**2 10 continue 20 continue ss=mbrun**2/2d0*lambda(s,dmW**2,mHp**2)*ss C...t-channel term tt=1d0/(t-mt**2)**2*( $ mt**2*mtrun**2/dtan(beta)**2*(2d0*dmW**2+pT**2)+ $ mbrun**2*dtan(beta)**2*(2d0*dmW**2*pT**2+t**2)) C...interferens term ts=0d0 do 30 i=1,3 ts=ts+dreal(ghichw(i)*ghibb(i)*prop(s,hkf(i)))/dcos(beta) 30 continue ts=mbrun**2*dtan(beta)/(t-mt**2)*(dmW**2*mHp**2-s*pT**2-t**2)*ts me_bbWmHp=Gf**2/(24d0*pi*s)*(ss+ts+tt) return end double precision function me_bbWpHm(s,t,pT,dmW,mHm) C...Double precision and integer declarations. IMPLICIT DOUBLE PRECISION(A-H, O-Z) IMPLICIT INTEGER(I-N) C...Parameters. COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) C...Particle properties + some flavour parameters. COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4) C...Supersymmetry parameters. COMMON/PYMSSM/IMSS(0:99),RMSS(0:99) double precision s,t,pT,lambda,beta,ss,ts double precision mbrun,mHm,dmW,mt,mtrun,Gf,pi,Q2 complex*16 ghichw,ghibb common/pybbwhg/ghichw(3),ghibb(3) complex*16 prop integer hkf(1:3) data hkf/25,35,36/ C...Propagator prop(s,kf)=1d0/dcmplx( $ s-pmas(kf,1)**2, $ pmas(kf,1)*pmas(kf,2)) lambda(x,y,z)=x**2+y**2+z**2-2*(x*y+y*z+z*x) beta=atan(RMSS(5)) Q2=(dmW+mHm)**2 mbrun=pymrun(5,Q2) mt=pmas(6,1) mtrun=pymrun(6,Q2) Gf=paru(1)*pyalem(Q2)/(sqrt(2d0)*paru(102)*pmas(24,1)**2) pi=paru(1) C...s-channel term ss=0 do 20 i=1,3 do 10 j=1,3 ss=ss+dreal( $ dconjg(ghichw(i))*ghichw(j)* $ prop(s,hkf(i))*dconjg(prop(s,hkf(j)))* $ dreal(conjg(ghibb(i))*ghibb(j)))/cos(beta)**2 10 continue 20 continue ss=mbrun**2/2d0*lambda(s,dmW**2,mHm**2)*ss C...t-channel term tt=1d0/(t-mt**2)**2*( $ mt**2*mtrun**2/dtan(beta)**2*(2d0*dmW**2+pT**2)+ $ mbrun**2*dtan(beta)**2*(2d0*dmW**2*pT**2+t**2)) C...interferens term ts=0 do 30 i=1,3 ts=ts+dreal(dconjg(ghichw(i))*dconjg(ghibb(i))*prop(s,hkf(i))) $ /dcos(beta) 30 continue ts=mbrun**2*dtan(beta)/(t-mt**2)*(dmW**2*mHm**2-s*pT**2-t**2)*ts me_bbWpHm=Gf**2/(24d0*pi*s)*(ss+ts+tt) return end