c======================================================================= c If you copy this code, please send a message to c RIDOLFI@VXCERN.CERN.CH c or to c MLM@VXCERN.CERN.CH c so that we can keep you informed of upgrades or modifications. c======================================================================= implicit real*8 (a-h,o-z) character *4 proc external radsig parameter(rmz=91.18d0,rmw=80.2d0,gz=2.49d0) common/beams/s,sthr common/prcpar/rma,rmh,rmhh,sbasq,proc common/bveg1/xl(10),xu(10),acc,ndim,ncall,itmx,nprn common/susy/rmt,rmb,rmq2,rmu2,rmd2,at,ab,rmu c Initialize parameters for Vegas do j=1,10 xl(j) = 0 xu(j) = 1 enddo acc = -1 ndim = 1 ncall = 10000 itmx = 10 1 write(*,*)'Enter sqrt(s), m(A), tan(beta), m(squark), proc' write(*,*)' proc = ''hz '' for h+Z production' write(*,*)' proc = ''fhz '' for h+ffbar production' write(*,*)' proc = ''hhz '' for H+Z production' write(*,*)' proc = ''fhhz'' for H+ffbar production' write(*,*)' proc = ''ha '' for h+A production' write(*,*)' proc = ''hha '' for H+A production' read(*,*)sr,rma,tb,rmsq,proc if(sr.lt.0)stop s=sr**2 rmt=176.d0 rmb=3.d0 q=80.d0 at=0.d0 ab=0.d0 rmu=0.d0 rmq2=rmsq**2 rmu2=rmsq**2 rmd2=rmsq**2 call susyhiggs(rma,tb,q,rmh,rmhh,rmhc,c2a,iflag) if (proc.eq.'hz ')then sthr=(rmh+rmz)**2 elseif(proc.eq.'fhz ')then sthr=rmh**2 elseif(proc.eq.'hhz ')then sthr=(rmhh+rmz)**2 elseif(proc.eq.'fhhz')then sthr=rmhh**2 elseif(proc.eq.'ha ')then sthr=(rmh+rma)**2 elseif(proc.eq.'hha ')then sthr=(rmhh+rma)**2 endif if(sthr.gt.s)then write(*,*)' unaccessible channel' goto 1 endif sb=tb/sqrt(1.d0+tb**2) cb=1.d0/sqrt(1.d0+tb**2) sa=-sqrt((1.d0-c2a)/2.d0) ca=sqrt((1.d0+c2a)/2.d0) sba=sb*ca-sa*cb sbasq=sba**2 c no isr xsect0=sig(s) c with isr call vegas(radsig,xsect1,sd,chisq) c output write(*,101) tb,rmsq write(*,102) rma,rmh,rmhh,rmhc write(*,*)proc,' xsec no isr =',xsect0,'pb' write(*,*)proc,' xsec with isr =',xsect1,'pb' if(sr.ne.srold) then write(1,100) sr srold=sr endif par=(1+rmsq+rmsq**2)*(tb+tb**2) if(par.ne.parold) then write(1,101) tb,rmsq parold=par endif write(1,102)rma,rmh,rmhh,rmhc write(1,103)proc,xsect0 write(1,104)proc,xsect1 goto 1 100 format(' Ecm = ',g10.4,'GeV') 101 format(' tan-beta=',g10.4,' m(squark)=',g10.4) 102 format(' m(A)=',g10.4,'GeV',' m(h)=',g10.4,'GeV', # ' m(H)=',g10.4,'GeV',' m(H+)=',g10.4,'GeV') 103 format(a,' xsec no isr (pb) =',g10.4) 104 format(a,' xsec with isr (pb) =',g10.4) goto 1 end function sig(s) implicit real*8 (a-z) character*4 proc external fh parameter(rmz=91.18d0,rmw=80.2d0,gz=2.49d0) common/prcpar/rma,rmh,rmhh,sbasq,proc common/beams/see,sthr common/mass/rm,ss data zero,uno,err/0.d0,1.d0,1.d-6/ ss=s aem=alfa(s) gf=1.16639d-5 hc2=0.3893857d9 pi=dacos(-1.d0) s2w=1-rmw**2/rmz**2 c2w=1-s2w clep=(0.5d0-s2w)**2+s2w**2 prop2=rmz**2/((s-rmz**2)**2+rmz**2*gz**2) coeff2=4*gf**2*rmz**3*gz/(3*pi**2) prop=1.d0/((s-rmz**2)**2+rmz**2*gz**2)/s**2 coeff=gf**2*rmz**4/(12*pi) if(s.lt.sthr)then sig=0 else if(proc.eq.'hz '.or.proc.eq.'fhz '.or.proc.eq.'ha ')then rm=rmh elseif(proc.eq.'hhz '.or.proc.eq.'fhhz'.or.proc.eq.'hha ')then rm=rmhh sbasq=1-sbasq endif c------------------------------------------------------- if(proc.eq.'hz '.or.proc.eq.'hhz ')then zlam=sqrt(s**2+rmz**4+rm**4-2*(s*rmz**2+rmz**2*rm**2+s*rm**2)) c eq. 4 bartl et al plb 237(90)229 sig=coeff*clep*prop*zlam*(12*rmz**2*s+zlam**2)*sbasq c------------------------------------------------------- elseif(proc.eq.'fhz '.or.proc.eq.'fhhz')then c eq. 5 bartl et al plb 237(90)229 sig=coeff2*clep*prop2*sbasq*dgauss(fh,zero,uno,err) c------------------------------------------------------- elseif(proc.eq.'ha '.or.proc.eq.'hha ')then alam=sqrt(s**2-2*s*(rm**2+rma**2)+(rm**2-rma**2)**2) sig=coeff*clep*prop*(1-sbasq)*alam**3 c------------------------------------------------------- endif endif sig=sig*hc2 end function fh(x) implicit real * 8 (a-z) parameter(rmz=91.18d0,rmw=80.2d0,gz=2.49d0) common/mass/rm,s xjac=1 csi1=atan((rmz**2-(sqrt(s)-rm)**2)/(rmz*gz)) csi2=atan(rmz/gz) csi=csi1+x*(csi2-csi1) xjac=xjac*(csi2-csi1) e=(rmz*gz*tan(csi)+s+rm**2-rmz**2)/(2*sqrt(s)) xjac=xjac*rmz*gz*(1+tan(csi)**2)/(2*sqrt(s)) num=(e**2-6*sqrt(s)*e+3*s+2*rm**2)*sqrt(e**2-rm**2) den=(s+rm**2-rmz**2-2*sqrt(s)*e)**2+rmz**2*gz**2 fh=num/den*xjac return end subroutine susyhiggs(ma,tb,q,mh,mhh,mhc,c2a,iflag) c --------------------------------------------------------------------- c | inputs: ma (CP-odd higgs mass) | c | tb (v2/v1) | c | q (scale) | c | common: mt,mb,mq2,mu2,md2,at,ab,mu | c | outputs: mh,mhh (light and heavy CP-even higgs masses) | c | mhc (charged higgs mass) | c | c2a (cos 2 alpha) | c | on output: iflag = 0 -> no unphysical squared mass | c | iflag = 1 -> squared stop mass < 2000 GeV**2 | c | iflag = 2 -> squared sbottom mass < 2000 GeV**2 | c | iflag = 3 -> negative squared h mass | c | iflag = 4 -> negative squared H+ mass | c --------------------------------------------------------------------- implicit real*8 (a-z) Integer Iflag external f1,f2,ff,f,g,g1,corr,h parameter(mw=80.d0,mz=91.174d0,alpha=1.d0/128.d0) common/susy/mt,mb,mq2,mu2,md2,at,ab,mu common/bcd/b12,c12,d12 common/quark/mt2,mb2 common/scale/q2 q2=q**2 mt2=mt**2 mb2=mb**2 mw2=mw**2 mz2=mz**2 iflag=0 pi=dacos(-1.d0) c2w=mw2/mz2 s2w=1-c2w t2w=s2w/c2w g2=4*pi*alpha/s2w sb=tb/dsqrt(1+tb**2) cb=1.d0/dsqrt(1+tb**2) s2b=2*sb*cb c2b=cb**2-sb**2 radt=3*g2*mt**4/(8*pi**2*sb**2*mw**2) radb=3*g2*mb**4/(8*pi**2*cb**2*mw**2) t1=(mq2+mu2+2*mt2+mz2*c2b/2.d0 # +sqrt((mq2-mu2+(8*mw2-5*mz2)*c2b/6.d0)**2 # +4*mt2*(at+mu/tb)**2))/2 t2=(mq2+mu2+2*mt2+mz2*c2b/2.d0 # -sqrt((mq2-mu2+(8*mw2-5*mz2)*c2b/6.d0)**2 # +4*mt2*(at+mu/tb)**2))/2 b1=(mq2+md2+2*mb2-mz2*c2b/2.d0 # +sqrt((mq2-md2-(4*mw2-mz2)*c2b/6.d0)**2 # +4*mb2*(ab+mu*tb)**2))/2 b2=(mq2+md2+2*mb2-mz2*c2b/2.d0 # -sqrt((mq2-md2-(4*mw2-mz2)*c2b/6.d0)**2 # +4*mb2*(ab+mu*tb)**2))/2 if(t2.lt.2000.d0)then iflag=1 goto 33 endif if(b2.lt.2000.d0)then iflag=2 goto 33 endif even11= mz**2*s2b/tb + ma**2*s2b*tb # + radb*dlog(b1*b2/mb**4) # + radb*2*ab*(ab+mu*tb)*f1(b1,b2) # + radt*mu**2*(at+mu/tb)**2*f2(t1,t2) # + radb*ab**2 *(ab+mu*tb)**2*f2(b1,b2) even22= mz**2*s2b*tb + ma**2*s2b/tb # + radt*dlog(t1*t2/mt**4) # + radt*2*at*(at+mu/tb)*f1(t1,t2) # + radt*at**2*(at+mu/tb)**2*f2(t1,t2) # + radb*mu**2*(ab+mu*tb)**2*f2(b1,b2) even12= - mz**2*s2b - ma**2*s2b # + radt*mu*(at+mu/tb)*f1(t1,t2) # + radb*mu*(ab+mu*tb)*f1(b1,b2) # + radt*at*mu*(at+mu/tb)**2*f2(t1,t2) # + radb*ab*mu*(ab+mu*tb)**2*f2(b1,b2) even11=.5d0*even11 even12=.5d0*even12 even22=.5d0*even22 mhsq =.5d0*(even11+even22-dsqrt((even11-even22)**2+4*even12**2)) mhhsq=.5d0*(even11+even22+dsqrt((even11-even22)**2+4*even12**2)) if(mhsq.lt.0)then iflag=3 goto 33 else mh =dsqrt(mhsq) mhh=dsqrt(mhhsq) endif c2a=(even11-even22)/dsqrt((even11-even22)**2+4*even12**2) b12 = -((cb**3*mw**4-cb*mb2*mw2)*sb**3-ab*mb2*mu*sb**2+(2*cb*mb 1 2*mt2-cb**3*mt2*mw2)*sb-at*cb**2*mt2*mu)/ 2.d0 c12 = ((2*cb**4*mw**6-2*cb**2*mb2*mw**4)*sb**6-3*ab*cb*mb2*mu*mw* 1 *2*sb**5+((2*cb**4*mb2-2*cb**4*mt2)*mw**4-2*cb**6*mw**6)*sb**4+ 2 (3*at*cb**3*mt2+3*ab*cb**3*mb2)*mu*mw2*sb**3+2*cb**6*mt2*mw** 3 4*sb**2-3*at*cb**5*mt2*mu*mw2*sb)*t2w c12 = c12-3*ab*cb*mb2*mu*mw2*sb**5+((-6*cb**4*mu2-6*cb**4*mt2- 1 6*cb**4*md2-6*cb**4*mb2)*mw**4+(6*cb**2*mb2*mu2+6*cb**2*mb2*mu* 2 *2+6*cb**2*mb2*mt2+6*cb**2*mb2*md2+6*cb**2*mb2**2-6*ab**2*cb**2 3 *mb2)*mw2-6*mb2**2*mu**2)*sb**4+((3*at*cb**3*mt2+3*ab*cb**3*m 4 b2)*mu*mw2+6*ab*cb*mb2*mu*mu2+(12*ab*cb*mb2*mt2+6*ab*cb*mb2*m 5 q2)*mu)*sb**3+((6*cb**4*mt2*mu2+6*cb**4*mt2*mu**2+6*cb**4*mt2** 6 2+(6*cb**4*md2+6*cb**4*mb2-6*at**2*cb**4)*mt2)*mw2-6*cb**2*mb 7 2*mt2*mu2-12*cb**2*mb2*mt2*mu**2-12*cb**2*mb2*mt2**2+(-12*cb**2 8 *mb2*mq2-6*cb**2*mb2*md2-12*cb**2*mb2**2+(6*at**2+12*ab*at+6*ab 9 **2)*cb**2*mb2)*mt2)*sb**2+((6*at*cb**3*mq2+6*at*cb**3*md2+12*a : t*cb**3*mb2)*mt2*mu-3*at*cb**5*mt2*mu*mw2)*sb-6*cb**4*mt2**2* ; mu**2 c12 = -c12/( 12.d0*cb*sb) d12 = ((-12*cb**4*mu2-12*cb**4*mt2+24*cb**4*md2+24*cb**4*mb2)*mw* 1 *6+(12*cb**2*mb2*mu2-24*cb**2*mb2*mu**2+12*cb**2*mb2*mt2-24*cb* 2 *2*mb2*md2-24*cb**2*mb2**2+24*ab**2*cb**2*mb2)*mw**4+24*mb2**2* 3 mu**2*mw2)*sb**6+((12*cb**6*mu2+12*cb**6*mt2-24*cb**6*md2-24* 4 cb**6*mb2)*mw**6+((12*cb**4*mt2-12*cb**4*mb2)*mu2+(12*cb**4*mt2 5 +24*cb**4*mb2)*mu**2+12*cb**4*mt2**2+(-24*cb**4*md2-36*cb**4*mb 6 2-12*at**2*cb**4)*mt2+24*cb**4*mb2*md2+24*cb**4*mb2**2-24*ab**2 7 *cb**4*mb2)*mw**4+(-12*cb**2*mb2*mt2*mu2+(-12*cb**2*mb2*mt2-24* 8 cb**2*mb2**2)*mu**2-18*cb**2*mb2*mt2**2+(-12*cb**2*mb2*mq2+24*c 9 b**2*mb2*md2+18*cb**2*mb2**2+(12*at**2+12*ab*at-24*ab**2)*cb**2 : *mb2)*mt2)*mw2)*sb**4 d12 = d12+12*ab*cb*mb2*mu*mw**4*sb**7+((6*at*cb**3*mt2-24*ab*cb* 1 *3*mb2)*mu*mw**4+(6*ab*cb*mb2*mu*mu2+(-18*ab*cb*mb2*mt2-24*ab*c 2 b*mb2*mq2)*mu)*mw2)*sb**5+((12*ab*cb**5*mb2-12*at*cb**5*mt2)* 3 mu*mw**4+(((12*at*cb**3*mq2+6*at*cb**3*md2+(18*at+18*ab)*cb**3* 4 mb2)*mt2+24*ab*cb**3*mb2*mq2)*mu-6*ab*cb**3*mb2*mu*mu2)*mw2)* 5 sb**3+((-12*cb**6*mt2*mu2-12*cb**6*mt2*mu**2-12*cb**6*mt2**2+(2 6 4*cb**6*md2+24*cb**6*mb2+12*at**2*cb**6)*mt2)*mw**4+(12*cb**4*m 7 b2*mt2*mu2+(12*cb**4*mb2*mt2-12*cb**4*mt2**2)*mu**2+18*cb**4*mb 8 2*mt2**2+(12*cb**4*mb2*mq2-24*cb**4*mb2*md2-18*cb**4*mb2**2+(-1 9 2*at**2-12*ab*at+24*ab**2)*cb**4*mb2)*mt2)*mw2)*sb**2+(6*at*c : b**7*mt2*mu*mw**4+(-12*at*cb**5*mq2-6*at*cb**5*md2-18*at*cb**5* ; mb2)*mt2*mu*mw2)*sb+12*cb**6*mt2**2*mu**2*mw2 d12 = t2w*d12 d12 = d12+((8*cb**4*mw**8-8*cb**2*mb2*mw**6)*sb**8-4*ab*cb*mb2*m 1 u*mw**4*sb**7+(-16*cb**6*mw**8+(16*cb**4*mb2-8*cb**4*mt2)*mw**6 2 +7*cb**2*mb2*mt2*mw**4)*sb**6+(2*at*cb**3*mt2+8*ab*cb**3*mb2)*m 3 u*mw**4*sb**5+(8*cb**8*mw**8+(16*cb**6*mt2-8*cb**6*mb2)*mw**6-1 4 4*cb**4*mb2*mt2*mw**4)*sb**4+(-4*at*cb**5*mt2-4*ab*cb**5*mb2)*m 5 u*mw**4*sb**3+(7*cb**6*mb2*mt2*mw**4-8*cb**8*mt2*mw**6)*sb**2+2 6 *at*cb**7*mt2*mu*mw**4*sb)*t2w**2 d12 = d12+(-18*ab*cb*mb2*mu*mu2-18*ab*cb*mb2*mt2*mu)*mw2*sb**5 1 +(9*cb**6*mb2*mt2*mw**4+((36*cb**4*md2+36*cb**4*mb2)*mt2*mu2+(3 2 6*cb**4*md2+72*cb**4*mb2)*mt2*mu**2+(36*cb**4*md2+54*cb**4*mb2) 3 *mt2**2+(-36*at**2*cb**4*md2-18*cb**4*mb2**2+(36*ab*at-36*at**2 4 )*cb**4*mb2)*mt2)*mw2+(-36*cb**2*mb2*md2-36*cb**2*mb2**2+36*a 5 b**2*cb**2*mb2)*mt2*mu2-36*cb**2*mb2*mt2*mu**4+((-72*cb**2*mb2* 6 mq2-72*cb**2*mb2**2-72*ab*at*cb**2*mb2)*mt2-72*cb**2*mb2*mt2**2 7 )*mu**2+(-36*cb**2*mb2*mq2-36*cb**2*mb2*md2-72*cb**2*mb2**2+36* 8 ab**2*cb**2*mb2)*mt2**2+(-36*cb**2*mb2*mq2**2+(72*ab*at*cb**2*m 9 b2-36*cb**2*mb2**2)*mq2+36*at**2*cb**2*mb2*md2+36*at**2*cb**2*m : b2**2-36*ab**2*at**2*cb**2*mb2)*mt2)*sb**2 d12 = d12+9*cb**2*mb2*mt2*mw**4*sb**6+(((-36*cb**4*md2-36*cb**4* 1 mb2)*mu2+(-36*cb**4*md2-54*cb**4*mb2)*mt2)*mw**4+((36*cb**2*mb2 2 *mu**2+36*cb**2*mb2*md2+36*cb**2*mb2**2-36*ab**2*cb**2*mb2)*mu2 3 +72*cb**2*mb2*mt2*mu**2-18*cb**2*mb2*mt2**2+(36*cb**2*mb2*md2+5 4 4*cb**2*mb2**2+(36*ab*at-36*ab**2)*cb**2*mb2)*mt2)*mw2-36*mb2 5 **2*mu**2*mu2-36*mb2**2*mt2*mu**2)*sb**4+((18*ab*cb**3*mb2*mu*m 6 u2+(18*at*cb**3*md2+(90*at+90*ab)*cb**3*mb2)*mt2*mu)*mw2+(36* 7 ab*cb*mb2*mt2+36*ab*cb*mb2*mq2)*mu*mu2-36*at*cb*mb2*mt2*mu**3+( 8 36*ab*cb*mb2*mt2**2+(36*ab*cb*mb2*mq2-72*at*cb*mb2**2-36*ab*at* 9 *2*cb*mb2)*mt2)*mu)*sb**3+((-18*at*cb**5*md2-18*at*cb**5*mb2)*m : t2*mu*mw2-36*ab*cb**3*mb2*mt2*mu**3+(((36*at*cb**3*md2+36*at* ; cb**3*mb2)*mq2+36*at*cb**3*mb2*md2+36*at*cb**3*mb2**2-36*ab**2* < at*cb**3*mb2)*mt2-72*ab*cb**3*mb2*mt2**2)*mu)*sb+(-36*cb**4*md2 = -36*cb**4*mb2)*mt2**2*mu**2 d12 = d12/( 72.d0*cb*sb) cost=3*g2/(16*pi**2*mw2*cb*sb) dterms=3*g2/(64*pi**2)*(mt2/sb**2+mb2/cb**2-mw2)* # (h(mt2)-h(mb2))/(mt2-mb2) kt=mt2*at*mu/(2*sb**2) kb=mb2*ab*mu/(2*cb**2) mhcsq=mw2+ma**2 # +3*g2/(8*pi**2*s2b*mw2*sb**2*cb**2)*corr(t1,b1,t2,b2) # +cost*(kt*f(t1,t2)+kb*f(b1,b2))-dterms if(mhcsq.ge.0)then mhc=dsqrt(mhcsq) else iflag=4 goto 33 endif 33 return end function f1(sqm1,sqm2) implicit real * 8 (a-h,o-z) if(sqm1.eq.sqm2)then f1=1.d0/sqm2 else f1=dlog(sqm1/sqm2)/(sqm1-sqm2) endif return end function f2(sqm1,sqm2) implicit real * 8 (a-h,o-z) if(sqm1.eq.sqm2)then f2=-1.d0/(6.d0*sqm2**2) else f2=(2.d0-(sqm1+sqm2)/(sqm1-sqm2)*dlog(sqm1/sqm2))/(sqm1-sqm2)**2 endif return end function corr(x1,x2,x3,x4) implicit real * 8 (a-h,o-z) if(x1.eq.x2.and.x2.eq.x3.and.x3.eq.x4)then corr=0 elseif(x1.eq.x3.and.x2.eq.x4)then corr=(g1(x1)+g1(x2))/(x1-x2)**2 elseif(x1.eq.x3)then corr=g1(x1)/((x1-x2)*(x1-x4))+g(x2)/((x2-x1)**2*(x2-x4)) # +g(x4)/((x4-x1)**2*(x4-x2)) elseif(x2.eq.x4)then corr=g1(x2)/((x2-x1)*(x2-x3))+g(x1)/((x1-x2)**2*(x1-x3)) # +g(x3)/((x3-x2)**2*(x3-x1)) else corr=g(x1)/((x1-x2)*(x1-x3)*(x1-x4)) # +g(x2)/((x2-x1)*(x2-x3)*(x2-x4)) # +g(x3)/((x3-x1)*(x3-x2)*(x3-x4)) # +g(x4)/((x4-x1)*(x4-x2)*(x4-x3)) endif return end function g(x) implicit real*8 (a-z) common/bcd/b12,c12,d12 g=-(b12*x**2+c12*x+d12)*ff(x) return end function g1(x) implicit real*8 (a-z) common/bcd/b12,c12,d12 g1=-(b12*x**2+c12*x+d12)*f(x,x)-(2*b12*x+c12)*ff(x) return end function f(x1,x2) implicit real*8 (a-z) common/quark/mt2,mb2 if(x1.eq.x2)then f=mt2/(mt2-mb2)*(log(x1/mt2)+1)-mb2/(mt2-mb2)*(log(x1/mb2)+1) else f=(ff(x1)-ff(x2))/(x1-x2) endif return end function ff(x) implicit real*8 (a-z) common/quark/mt2,mb2 ff=x*(mt2*log(x/mt2)-mb2*log(x/mb2))/(mt2-mb2) return end function h(x) implicit real*8 (a-z) common/scale/q2 if(x.eq.0)then h=0 else h=2*x*(log(x/q2)-1) endif return end