real function contrainte2(x) vector eta vector rho vector eta_up vector babar vector eta_amb vector rho_amb common/fit_param/cenvbc,sigvbc,cenvubcb,sigvubcb,cendelmd &,sigdelmd,cenepsK,sigepsK,cenmt,sigmt,cenbk,sigbk, &cenfbb,sigfbb,cenxi_s,sigxi_s, &cens2a,sigs2a,cens2b,sigs2b common/cte_val/pi,h,gf,xmw,xmbd,xmbs,xmk,delmk,xmc,xlamb, &xnu1,xnu2,xnu3,xnub,fk logical babar_flag,dms_flag common /fit_drive/babar_flag,dms_flag etamax = eta_up(1) *- Initialisation : if (x.eq.0.) then pi = acos(-1.) h = 6.582122e-13 gf = 1.16639e-5 ! Fermi Constant xmw = 80.33 ! W Mass ( error neglected ) xmbd = 5.279 ! B0d Mass ( error neglected ) xmbs = 5.369 ! B0s Mass ( error neglected ) xmk = .4977 ! K0 Mass ( error neglected ) delmk = 3.510e-15 ! DmK (error neglected) xmc = 1.3 ! c quark Mass ( error neglected ) xlamb = .2205 ! Lambda ( error = .0018 neglected ) xnu1 = 1.38 ! eta parameters for EpsK (Buras) xnu2 = .57 xnu3 = .47 xnub = .55 fk = .161 ! K decay constant *-------------------------------------------------------------------------- open (UNIT=10,FILE='parameters.input',STATUS='OLD') c read(10,1010)flg_name,dms_flag read(10,1010)flg_name,babar_flag 1010 format(a14,1x,l1) read(10,*) c c Input values c read(10,1020)var_name,cenvbc,sigvbc read(10,1020)var_name,cenvubcb,sigvubcb read(10,1020)var_name,cendelmd,sigdelmd read(10,1020)var_name,cenepsK,sigepsK cenepsK = cenepsK*1.e-03 sigepsK = sigepsK*1.e-03 read(10,1020)var_name,cenmt,sigmt read(10,1020)var_name,cenbk,sigbk read(10,1020)var_name,cenfbb,sigfbb read(10,1020)var_name,cenxi_s,sigxi_s read(10,1020)var_name,cens2a,sigs2a read(10,1020)var_name,cens2b,sigs2b 1020 format(1x,a9,2(2x,f8.4)) close(10) if (babar_flag) babar(1) = 1. endif c c Dmd constraint c if(x.eq.1.) then xt = (cenmt/xmw)**2 C1 = gf**2/(6*pi**2)*xlamb**2*xmw**2*(xmbd/h)*xnub*si(xt) dmd = cendelmd +2.*sigdelmd vcb = cenvbc - 2.*sigvbc fbb = cenfbb - sigfbb C2max = dmd/(C1*vcb**2*fbb**2) dmd = cendelmd - 2.*sigdelmd vcb = cenvbc + 2.*sigvbc fbb = cenfbb + sigfbb C2min = dmd/(C1*vcb**2*fbb**2) xmin = 1.-sqrt(C2min) xmax = 1. step = abs(xmax-xmin)/50. do i = 1,50 x = xmin + float(i-1)*step rho(i) = x eta(i) = dmd_fc(x,c2min,etamax) enddo rho(51) = 1. eta(51) = dmd_fc(1.,c2min,etamax) ilast = 51 xmin = 1 xmax = 1. - sqrt(c2max) step = abs(xmax-xmin)/48. do i = 1,48 x = xmin - float(i-1)*step rho(i+ilast) = x eta(i+ilast) = dmd_fc(x,c2max,etamax) enddo x = 1. - sqrt(c2max) rho(100) = x eta(100) = dmd_fc(x,c2max,etamax) endif c c Epsilon K constraint c if(x.eq.2.) then epsK = cenepsK xt = (cenmt/xmw)**2 xc = (xmc/xmw)**2 C1 = gf**2*fk**2*xmk*xmw**2/(6.*sqrt(2.)*pi**2*delmk) B1_f0 = -xnu1*si(xc) + xnu3*sij(xc,xt) B2_f0 = vcb**2*xnu2*si(xt) vcb = cenvbc + 2.*sigvbc bk = cenbk + sigbk C2max = C1*BK*vcb**2*xlamb**2 B2_f0max = vcb**2*xnu2*si(xt) vcb = cenvbc - 2.*sigvbc bk = cenbk - sigbk C2min = C1*BK*vcb**2*xlamb**2 B2_f0min = vcb**2*xnu2*si(xt) xmin = -1. xmax = 1. step = abs(xmax-xmin)/50. do i = 1,50 x = xmin + float(i-1)*step y = fc_eps(epsk,C2max,B1_f0,b2_f0max,x,etamax) rho(i) = x eta(i) = y enddo x = 1. y = fc_eps(epsk,C2max,B1_f0,b2_f0max,x,etamax) rho(51) = x eta(51) = y y = fc_eps(epsk,C2min,B1_f0,b2_f0min,x,etamax) rho(52) = x eta(52) = y ilast = 52 xmin = 1 xmax = -1. step = abs(xmax-xmin)/48. do i = 1,48 x = xmin - float(i)*step y = fc_eps(epsk,C2min,B1_f0,b2_f0min,x,etamax) rho(ilast+i) = x eta(ilast+i) = y enddo endif c c Vub/Vcb constraint c if(x.eq.3.) then valm = cenvubcb - 2.*sigvubcb valp = cenvubcb + 2.*sigvubcb xmin = -valp/xlamb xmax = -valm/xlamb step = abs(valp/xlamb-valm/xlamb)/10. do i = 1,10 x = xmin + float(i)*step rho(i) = x eta(i) = 0. ilast = i enddo xmin = -valm/xlamb xmax = valm/xlamb step=abs(xmax-xmin)/30. do i = 1,30 x = xmin + float(i)*step rho(ilast+i)=x y = (valm/xlamb)**2-x**2 if (y.le.0) then y = 0. else y = sqrt(y) endif eta(ilast+i)=y enddo ilast = 40 xmin = valm/xlamb xmax = valp/xlamb step=abs(xmax-xmin)/10. do i=1,10 x = xmin + float(i)*step rho(ilast+i) = x eta(ilast+ i) = 0. enddo ilast = 50 xmin = -valp/xlamb xmax = valp/xlamb step = abs(xmax-xmin)/50. do i = 1,50 x = xmax- float(i)*step rho(ilast+i)=x y = (valp/xlamb)**2-x**2 if (y.le.0) then y = 0. else y = sqrt(y) endif eta(ilast+i)=y enddo endif c c Sin(2A) Sin(2B) constraints c if(x.eq.4.) then nx = 100 ny = 100 deltarho=2./nx deltaeta=1./ny *scan the (rho,eta) plane do irho=1,nx rhov =-1 -deltarho/2. +irho*deltarho do ieta=1,ny etav=-deltaeta/2. +ieta*deltaeta *sin(2alpha) ca1 = etav/sqrt(rhov**2+etav**2) sa1 = rhov/sqrt(rhov**2+etav**2) ca2 = etav/sqrt((1.-rhov)**2+etav**2) sa2 = (1.-rhov)/sqrt((1.-rhov)**2+etav**2) s12=(sa1*ca2 + ca1*sa2) c12=(ca1*ca2 - sa1*sa2) s2a = 2.*s12*c12 *sin(2beta) s2b=(2.*etav*(1-rhov))/((1-rhov)**2+etav**2) chi2_s2a = ((s2a-cens2a)/sigs2a)**2 chi2_s2b = ((s2b-cens2b)/sigs2b)**2 chiii=chi2_s2a + chi2_s2b call hfill(2122,rhov,etav,chiii) enddo enddo alpha0=.5*asin(cens2a) beta0= .5*asin(cens2b) call getrhoeta(alpha0,beta0,rho_amb(1),eta_amb(1)) alpha=pi/2.-alpha0 beta=beta0 call getrhoeta(alpha,beta,rho_amb(2),eta_amb(2)) alpha=alpha0 beta=pi/2.-beta0 call getrhoeta(alpha,beta,rho_amb(3),eta_amb(3)) alpha=pi/2.-alpha0 beta=pi/2.-beta0 call getrhoeta(alpha,beta,rho_amb(4),eta_amb(4)) endif 999 continue return end ************************************************************* SUBROUTINE GETRHOETA(ALPHA,BETA,RHO,ETA) rho=tan(beta)/(tan(beta)-tan(alpha+beta)) eta=-rho*tan(alpha+beta) RETURN END ************************************************************* real function rplus(z) rplus = (1./z)*(1.+sqrt(1-z**2)) return end ************************************************************* real function rminus(z) rminus = (1./z)*(1.-sqrt(1-z**2)) return end ************************************************************* function fc_eps(epsk,C2,B1_f0,b2_f0,x,etamax) y = epsk/(C2*(B1_f0+B2_f0*(1.-x))) if (y.gt.etamax) y = etamax if( y .lt.0.) y = 0. fc_eps = y return end ************************************************************* function dmd_fc(x,cte,etamax) y = cte - (1.-x)**2 if (y.le.0) then y = 0. else y = sqrt(y) endif if(y.gt.etamax) y = etamax dmd_fc = y return end ************************************************************* FUNCTION si(x) si = 1./4.+9./4.*1./(1.-x)-3./2.*1./(1.-x)**2 si = si-3./2.*x**2/(1.-x)**3*alog(x) si = x*si return end ************************************************************* FUNCTION sij(xi,xj) sij1 = (1./4.+3./2.*1./(1.-xi)-3./4.*1./(1.-xi)**2)* & alog(xi)/(xi-xj) sij2 = (1./4.+3./2.*1./(1.-xj)-3./4.*1./(1.-xj)**2)* & alog(xj)/(xj-xi) sij3 = -3./4.*1./((1.-xi)*(1.-xj)) sij = xi*xj*(sij1+sij2+sij3) return end