c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: subroutine dilep(p_ele,b0ele,p_muo,b0muo,output) c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c subroutine calculates probability of a hypothesis that a lepton and c a jet have originated from t->Wb; W->lepton+neutrino, b->jet cascade c decay together with another t-quark decaying via t->Wb; W->anti-lepton+nu, c b->jet channel. Probability spanning the multidimensional parameter space c is projected on the top quark mass axis, thus representing the result c as the probability of a hypothesis as a function of the top quark mass. c Mass of the top and anti-top quarks are required to be the same. c c include files c ------------- include 'hadtop_fast.inc' include 'x1x2.inc' include 'ellipse_e.inc' include 'ellipse_m.inc' include 'partons_dilep.inc' include 'qq_talk.inc' c local declarations c ------------------ integer i,j,k, ii,jj,kk, im,ib,ibb integer ibin integer nofit, ifault real p_ele(4), p_muo(4),jet_ib(4), neu_te(3),neu_tm(3) real output(200) real xmt, pt_ttb,prob_ttb real x1,x2, qq1,qq2,gg1,gg2, pgg,pqq, xprob real b0ele(4),b0muo(4),theta,pjet integer iter real b_ele(4),b_ele_corr(4),sig_be,pt_be,eps_be real b_muo(4),b_muo_corr(4),sig_bm,pt_bm,eps_bm real probl,probe,probm, prob_be,prob_bm, prob_all, prob_cur logical first real elsigma,musigma C integer bj_grid,bc_grid,bn_grid C integer wm_grid,wc_grid,wn_grid C integer el_grid,ec_grid,en_grid c Lorentz boosts arrays c --------------------- real lep_4a(4),lep_4b(4),top_4a(4),top_4b(4) integer iboost real probl_cms data first /.true./ data cone_size /1.0/ c c PDFLIB c ------ real*8 rpdf character cpdf*20 real*8 dtop_m,dx1,dx2 real*8 upv1,dnv1,sea1,str1,chm1,bot1,top1,gl1 real*8 upv2,dnv2,sea2,str2,chm2,bot2,top2,gl2 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C No More Declarations CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C bc_grid=(bj_grid+1)/2 C bn_grid=(bj_grid-1)/2 C wc_grid=(wm_grid+1)/2 C wn_grid=(wm_grid-1)/2 C ec_grid=(el_grid+1)/2 C en_grid=(el_grid-1)/2 c++ 3-feb-2011 write(63,*) 'P_ele: ',p_ele write(63,*) 'b0ele: ',b0ele write(63,*) 'P_muo: ',p_muo write(63,*) 'b0muo: ',b0muo c++ 3-feb-2011 c PDFLIB initialization c --------------------- if(first) then cpdf = 'NPTYPE' !Particle type. rpdf = 1 !1=nucleon. CALL PDFSET(CPDF,RPDF) c++ 4-feb-2011 write(63,*) 'P_ele: ',p_ele write(63,*) 'b0ele: ',b0ele write(63,*) 'P_muo: ',p_muo write(63,*) 'b0muo: ',b0muo c++ 4-feb-2011 c c cpdf = 'NGROUP' !Author group. c rpdf = 1 !1=DO. c CALL PDFSET(CPDF,RPDF) c cpdf = 'NSET' !pdf set. c rpdf = 7 !7=DO2 c CALL PDFSET(CPDF,RPDF) c c MRS parametrization c cpdf = 'NGROUP' !Author group. rpdf = 3 !3=MRS. CALL PDFSET(CPDF,RPDF) cpdf = 'NSET' !pdf set. rpdf = 30 !30=MRS D0' (MS_bar) c rpdf = 33 !33=MRS D0' (DIS) c rpdf = 39 !39=MRS (Ap). c rpdf = 40 !40=MRS (Ap) Fit. CALL PDFSET(CPDF,RPDF) first=.false. endif c reset the array into which the probability will be projected c (mt_bin bins from mt_low to mt_low+mt_bin*mt_gap GeV) do i=1,mt_bin sum(i)=0.0 enddo c reset other counters x_1=0.0 x_2=0.0 x12_p=0.0 p_hadt=0.0 p_lepb=0.0 p_ttb=0.0 topmin_e=0.0 topmin_m=0.0 atmin_e=1.e3 atmin_m=1.e3 prob_all=0.0 c prepare to store the partons 4-momenta neu_e(1)=0.0 neu_e(2)=0.0 neu_e(3)=0.0 neu_m(1)=0.0 neu_m(2)=0.0 neu_m(3)=0.0 bj_ele(1)=0.0 bj_ele(2)=0.0 bj_ele(3)=0.0 bj_ele(4)=0.0 bj_muo(1)=0.0 bj_muo(2)=0.0 bj_muo(3)=0.0 bj_muo(4)=0.0 e_lep(1)=0.0 e_lep(2)=0.0 e_lep(3)=0.0 m_lep(1)=0.0 m_lep(2)=0.0 m_lep(3)=0.0 weights=0.0 c index will count the number of successful pairs for ptenu+b momentum space c probability of lepton and b-quark jet candidate to originate from t decay c (Dalitz+Goldstein Phys.Rev.D) call ellipse_e(p_ele,b_ele_corr,nofit) if(nofit.eq.1) go to 1300 c check how this top quark matches the other top quark in the transverse c momentum plane (in this case both deacy leptonically) do 1200 k=1,el_grid ! ellipse for t->Wb, W->enu top top_e(1)=top_ele(k,1) top_e(2)=top_ele(k,2) top_e(3)=top_ele(k,3) top_e(4)=sqrt(top_m**2+top_e(1)**2+top_e(2)**2 < +top_e(3)**2) c neutrino neu_te(1)=top_e(1)-p_ele(1)-b_ele_corr(1) neu_te(2)=top_e(2)-p_ele(2)-b_ele_corr(2) neu_te(3)=top_e(3)-p_ele(3)-b_ele_corr(3) c "wiggle" the b_muo jet do 2300 ii=1,bj_grid ! b-jet accompanying the W->munu decay eps_bm=float(ii-bc_grid)/float(bn_grid) eps_bm=eps_bm*nsigma*sig_bm b_muo_corr(4)=b_muo(4)+eps_bm c keep energy positive for large errors c write(63,*) b_muo_corr(4) if(b_muo_corr(4).le.0) go to 2300 c find (Gaussian) probability of such a deviation from the initial b-jet call gauxx(b_muo_corr(4),b_muo(4),sig_bm,prob_bm) eps_bm=b_muo_corr(4)/b_muo(4) do jj=1,3 b_muo_corr(jj)=b_muo(jj)*eps_bm enddo ! end jj c write(63,*) ' b_mu: ', b_muo c write(63,*) ' corr: ', b_muo_corr c find the ellipse in the t->munu+b momentum space c probability of lepton and b-quark jet candidate to originate from t decay c (Dalitz+Goldstein Phys.Rev.D) call ellipse_m(p_muo,b_muo_corr,nofit) if(nofit.eq.1) go to 2300 c check how this top quark matches the other top quark in the transverse c momentum plane (in this case both deacy leptonically) do 2200 kk=1,el_grid ! ellipse for t->Wb, W->munu top top_mu(1)=top_muo(kk,1) top_mu(2)=top_muo(kk,2) top_mu(3)=top_muo(kk,3) top_mu(4)=sqrt(top_m**2+top_mu(1)**2+top_mu(2)**2 < +top_mu(3)**2) c neutrino neu_tm(1)=top_mu(1)-p_muo(1)-b_muo_corr(1) neu_tm(2)=top_mu(2)-p_muo(2)-b_muo_corr(2) neu_tm(3)=top_mu(3)-p_muo(3)-b_muo_corr(3) c find the ttbar transverse momentum vector pt_ttb=sqrt((top_e(1)+top_mu(1))**2+(top_e(2) < +top_mu(2))**2) c write(63,*) 'top_e: ',top_e,top_m c write(63,*) 'top_m: ',top_mu,pt_ttb c ISAJET ttbar pt distribution (Kuni Kondo's fit) c K Sliwa's fit to VECBOS xmt=pt_ttb/top_m if(xmt.gt.fx_cut) go to 2200 cKuni prob_ttb=128.95*xmt*exp(6.7*xmt**2-13.39*xmt) prob_ttb=611.*xmt*exp(21.*xmt**2-28.*xmt) c find the probability for this configuration of lepton+4jets to have c originated from t+tbar production with x1 and x2 postulated here x1=(top_mu(4)+top_e(4)+top_mu(3)+top_e(3))/TEV_Energy x2=(top_mu(4)+top_e(4)-top_mu(3)-top_e(3))/TEV_Energy c Insert the Duke-Owens functions here, as in Dalitz+Goldstein Phys.Rev.D c article (borrowed from the code written in the spring 1992 by Gary and I) if(x1.lt.0.01) x1=0.01 if(x2.lt.0.01) x2=0.01 if(x1.ge.1.00) then qq1=0.0 gg1=0.0 go to 2200 endif if(x2.ge.1.00) then qq2=0.0 gg2=0.0 go to 2200 endif c Use both q+qbar and gluon fusion into t+tbar to get probability pqq=0.0 pgg=0.0 call cross_dilep(x1,x2,pqq,pgg) c use PDFLIB structure function parametrizations dtop_m=top_m dx1=x1 dx2=x2 call structf(dx1,dtop_m,upv1,dnv1,sea1,str1, < chm1,bot1,top1,gl1) call structf(dx2,dtop_m,upv2,dnv2,sea2,str2, < chm2,bot2,top2,gl2) qq1=(upv1+dnv1+sea1+str1+chm1+bot1)/dx1 qq2=(upv2+dnv2+sea2+str2+chm2+bot2)/dx2 gg1=gl1/dx1 gg2=gl2/dx2 c write(65,*) 'PDF ',top_m,x1,x2,qq1,qq2,gg1,gg2 xprob=pqq*qq1*qq2+pgg*gg1*gg2 c c call hf2(7001,top_m,pqq,1.0) c call hf2(7002,top_m,pgg,1.0) c call hf2(7101,x1,qq1,1.0) c call hf2(7102,x2,qq2,1.0) c call hf2(7201,x1,gg1,1.0) c call hf2(7202,x2,gg2,1.0) c call hf2(7300,top_m,xprob,1.0) c W->enu decay c boost lepton to t cms do iboost=1,4 top_4a(iboost)=-top_e(iboost) enddo top_4a(4)=top_e(4) c lepton do iboost=1,4 lep_4a(iboost)=p_ele(iboost) lep_4b(iboost)=0.0 enddo call lorenb(top_m,top_4a,lep_4a,lep_4b) c calculate probl factor in t cms fact=lep_4b(4)/top_m probl_cms=24.*fact*(1.-2.*fact) probe=probl_cms c probe=4*(bdotl_e+w_mass**2/2.)* c < (top_m**2-b_mass**2-2.*bdotl_e-w_mass**2)/ c < ((top_m**2-b_mass**2)**2+ c < (top_m**2+b_mass**2)*(w_mass**2)-2*w_mass**4) if(probe.lt.0) go to 2200 c W->munu decay c boost lepton to t cms do iboost=1,4 top_4a(iboost)=-top_mu(iboost) enddo top_4a(4)=top_mu(4) c lepton do iboost=1,4 lep_4a(iboost)=p_muo(iboost) lep_4b(iboost)=0.0 enddo call lorenb(top_m,top_4a,lep_4a,lep_4b) c calculate probl factor in t cms fact=lep_4b(4)/top_m probl_cms=24.*fact*(1.-2.*fact) probm=probl_cms c probm=4*(bdotl_m+w_mass**2/2.)* c < (top_m**2-b_mass**2-2.*bdotl_m-w_mass**2)/ c < ((top_m**2-b_mass**2)**2+ c < (top_m**2+b_mass**2)*(w_mass**2)-2*w_mass**4) if(probm.lt.0) go to 2200 c c call hf1(7310,probe,1.0) c call hf1(7311,probm,1.0) c call hf1(7312,prob_be,1.0) c call hf1(7313,prob_bm,1.0) prob_cur=prob_ttb*prob_be*prob_bm*probe*probm*xprob p_lepb=p_lepb+prob_cur*probm p_hadt=p_hadt+prob_cur*probe x_1=x_1+prob_cur*x1 x_2=x_2+prob_cur*x2 p_ttb=p_ttb+prob_cur*prob_ttb x12_p=x12_p+prob_cur*xprob c neutrino/weighted weights=weights+prob_cur neu_e(1)=neu_te(1)*prob_cur+neu_e(1) neu_e(2)=neu_te(2)*prob_cur+neu_e(2) neu_e(3)=neu_te(3)*prob_cur+neu_e(3) neu_m(1)=neu_tm(1)*prob_cur+neu_m(1) neu_m(2)=neu_tm(2)*prob_cur+neu_m(2) neu_m(3)=neu_tm(3)*prob_cur+neu_m(3) c b from leptonic decay bj_ele(1)=b_ele_corr(1)*prob_cur+bj_ele(1) bj_ele(2)=b_ele_corr(2)*prob_cur+bj_ele(2) bj_ele(3)=b_ele_corr(3)*prob_cur+bj_ele(3) bj_ele(4)=b_ele_corr(4)*prob_cur+bj_ele(4) c b from leptonic decay bj_muo(1)=b_muo_corr(1)*prob_cur+bj_muo(1) bj_muo(2)=b_muo_corr(2)*prob_cur+bj_muo(2) bj_muo(3)=b_muo_corr(3)*prob_cur+bj_muo(3) bj_muo(4)=b_muo_corr(4)*prob_cur+bj_muo(4) c electron e_lep(1)=p_ele(1)*prob_cur+e_lep(1) e_lep(2)=p_ele(2)*prob_cur+e_lep(2) e_lep(3)=p_ele(3)*prob_cur+e_lep(3) c muon m_lep(1)=p_muo(1)*prob_cur+m_lep(1) m_lep(2)=p_muo(2)*prob_cur+m_lep(2) m_lep(3)=p_muo(3)*prob_cur+m_lep(3) c lower bounds for Mt from lepton+jets topmin_e=tmin_e*prob_cur+topmin_e topmin_m=tmin_m*prob_cur+topmin_m atmin_e=amin1(atmin_e,tmin_e) atmin_m=amin1(atmin_m,tmin_m) c total prob_all=prob_cur+prob_all index=index+1 2200 continue ! end of loop over muon ellipse 2300 continue ! end of muon b-jet smearing loop 1200 continue ! loop over all points on the electron ellipse 1300 continue ! b-quark jet from the electron t decay c store the summed probabilities for this combination of lepton+4jets ibin=(top_m-mt_low)/mt_gap+1 if(ibin.lt.1) ibin=1 if(ibin.gt.mt_bin) ibin=mt_bin sum(ibin)=sum(ibin)+prob_all 1333 continue enddo ! Mt loop if(weights.gt.0) then c write(63,900) atmin_e,atmin_m,sum CC write(63,900) atmin_e,atmin_m endif 900 format(108(1h.),/,'minimum M(top) [GeV]/1/2 :',2f10.1) c 900 format(108(1h.),/,'minimum M(top) [GeV]/1/2 :',2f10.1,/,(10g10.2) output(1)=atmin_e output(2)=atmin_m output(3)=topmin_e output(4)=topmin_m output(5)=weights output(6)=index output(7)=x_1 output(8)=x_2 output(9)=x12_p output(10)=p_lepb output(11)=p_hadt output(12)=p_ttb do iter=1,3 output(12+iter)=e_lep(iter) output(16+iter)=m_lep(iter) output(20+iter)=neu_e(iter) output(24+iter)=neu_m(iter) output(28+iter)=bj_ele(iter) output(32+iter)=bj_muo(iter) enddo output(16)=sqrt(e_lep(1)**2 + e_lep(2)**2 + e_lep(3)**2) output(20)=sqrt(m_lep(1)**2 + m_lep(2)**2 + m_lep(3)**2) output(24)=sqrt(neu_e(1)**2 + neu_e(2)**2 + neu_e(3)**2) output(28)=sqrt(neu_m(1)**2 + neu_m(2)**2 + neu_m(3)**2) output(32)=bj_ele(4) output(36)=bj_muo(4) output(37)=cmfrnc_e output(38)=cmfrnc_m output(39)=elsigma output(40)=musigma output(41)=mt_bin do iter=1,mt_bin output(41+iter)=sum(iter) enddo return end SUBROUTINE MSIGMA(E,SIGMA) c Returns error on jet, based on E and Eta c Note Eta not implemented yet c Input E, energy and ETA detector eta of the jet REAL E,SIGMA c approximate statement functions REAL SIGT,SIG,SIGL DATA SIG/0.10/ SIGT(E) = 0.29856-0.68355E-02*E+0.85703E-04*E*E & -0.36847E-06*E*E*E SIGL(E) = 0.85/SQRT(E) IF(E.GT.100.) THEN SIGMA = SIG*E ELSE IF(E.LT.20) THEN SIGMA = SIGL(E)*E ELSE SIGMA =SIGT(E)*E ENDIF RETURN END subroutine gauxx(en,ej,sigma,pgauss) c Normalized Gaussian function real en,ej,sigma,pi,pgauss data pi /3.14159265359/ c pgauss=(exp(-(((en-ej)/sigma)**2)/2))/(sigma*sqrt(2.0*pi)) return end subroutine ellipse_e(p_lep,b_ele,nofit) c include files c ------------- include 'hadtop_fast.inc' include 'ellipse_e.inc' include 'x1x2.inc' c local declarations c ------------------ integer i,j,k real mtmin1 real m0s,b1,aa,kos,lam,zmt,al2(4),al2s,al1(4),alp,xp,yp,zp real ts,pi,circum_e real p_lep(4),b_ele(4) integer nofit data pi /3.14159265359/ c nofit=0 c p_be=0.0 p_e=0.0 do j=1,3 p_e=p_e+p_lep(j)**2 p_be=p_be+b_ele(j)**2 enddo ! end j en_e=p_lep(4) c use the true b-jet energy en_be=b_ele(4) ! not sure how important c en_be=sqrt(p_be+b_mass**2) ! difference it may cause bdotl_e=en_be*en_e bscal_e=0.0 do j=1,3 bscal_e=bscal_e+p_lep(j)*b_ele(j) bdotl_e=bdotl_e-p_lep(j)*b_ele(j) enddo ! end j c write(63,*) 'b.l; bscal_e',bdotl_e,bscal_e 60 mtmin1=sqrt((b_mass**2+2.*bdotl_e)*(w_mass**2+2.*bdotl_e)/ < (2.*bdotl_e)) tmin_e=mtmin1 c write(63,*) 'mt_min, top_e: ',mtmin1,top_m 61 if(mtmin1.ge.top_m) then c write(63,*) 'Mtmin exceeds top_m' nofit=1 goto 9009 endif 70 m0s=b_mass**2+2.*bdotl_e+(w_mass**2)*(1.+(en_be/en_e)- < bdotl_e/(2.*en_e**2)) c 100 write(63,*) 'm0=',sqrt(m0s),'mtmin=',mtmin1 200 b1=sqrt(p_be-(bscal_e/en_e)**2) 210 aa=bdotl_e/en_e c 220 write(63,*) 'ellipse plane angle = ',atan(b1/aa)*180./pi 230 kos=aa/sqrt(aa**2+b1**2) 300 lam=w_mass*sqrt((top_m**2-mtmin1**2)/(2.*bdotl_e)) c write(63,*) 'top_m=',top_m,'major _',lam/kos,'minor _',lam circum_e=2*pi*lam/sqrt(kos) cmfrnc_e=circum_e 500 zmt=(top_m**2-m0s)*kos/(2*aa) 1100 al2(1)=b_ele(2)*p_lep(3)-b_ele(3)*p_lep(2) 1101 al2(2)=b_ele(3)*p_lep(1)-b_ele(1)*p_lep(3) 1102 al2(3)=b_ele(1)*p_lep(2)-b_ele(2)*p_lep(1) 1103 al2s=0 1110 do j=1,3 1111 al2s=al2(j)**2+al2s 1112 enddo ! end j c write(63,*)'al2 & al2s',al2(1),al2(2),al2(3),al2s 1120 do j=1,3 1121 al2(j)=al2(j)/sqrt(al2s) 1122 enddo ! end j 1130 al1(1)=(al2(2)*p_lep(3)-al2(3)*p_lep(2))/en_e 1131 al1(2)=(al2(3)*p_lep(1)-al2(1)*p_lep(3))/en_e 1132 al1(3)=(al2(1)*p_lep(2)-al2(2)*p_lep(1))/en_e c write(63,*) 'eta_e','t(x)','t(y)','t(z)' 1200 do k=1,el_grid 1210 alp=(k-1)*2.*pi/(el_grid) eta_e(k)=(k-1)*360./(el_grid) 1220 xp=lam*cos(alp)+(w_mass**2)*(b1/aa)/(2*en_e) 1225 yp=lam*sin(alp) 1230 zp=xp*(b1/aa)+zmt/kos 1300 do j=1,3 1315 top_ele(k,j)=b_ele(j)+(zp+en_e-(w_mass**2/(4.*en_e)))* < p_lep(j)/en_e-xp*al1(j)-yp*al2(j) 1320 enddo ! end j 1390 ts=0 1400 do j=1,3 1410 ts=top_ele(k,j)**2+ts 1420 enddo ! end j c write(63,*) eta_e(k),top_ele(k,1),top_ele(k,2),top_ele(k,3) 2000 enddo ! end k 9009 return end subroutine ellipse_m(p_muo,b_muo,nofit) c include files c ------------- include 'hadtop_fast.inc' include 'ellipse_m.inc' include 'x1x2.inc' c local declarations c ------------------ integer i,j,k real mtmin1 real m0s,b1,aa,kos,lam,zmt,al2(4),al2s,al1(4),alp,xp,yp,zp real ts,pi,circum_m real p_muo(4),b_muo(4) integer nofit data pi /3.14159265359/ c nofit=0 c c write(63,*) 'p_muo: ', p_muo c write(63,*) 'b_muo: ', b_muo c p_bm=0.0 p_m=0.0 do j=1,3 p_m=p_m+p_muo(j)**2 p_bm=p_bm+b_muo(j)**2 enddo ! end j en_m=p_muo(4) c use the true b-jet energy en_bm=b_muo(4) ! not sure how important c en_bm=sqrt(p_bm+b_mass**2) ! difference it may cause bdotl_m=en_bm*en_m bscal_m=0.0 do j=1,3 bscal_m=bscal_m+p_muo(j)*b_muo(j) bdotl_m=bdotl_m-p_muo(j)*b_muo(j) enddo ! end j c write(63,*) 'b.l; bscal_m',bdotl_m,bscal_m 60 mtmin1=sqrt((b_mass**2+2.*bdotl_m)*(w_mass**2+2.*bdotl_m)/ < (2.*bdotl_m)) tmin_m=mtmin1 c write(63,*) 'mt_min, top_m: ',mtmin1,top_m 61 if(mtmin1.ge.top_m) then c write(63,*) 'Mtmin exceeds top_m' nofit=1 goto 9009 endif 70 m0s=b_mass**2+2.*bdotl_m+(w_mass**2)*(1.+(en_bm/en_m)- < bdotl_m/(2.*en_m**2)) c 100 write(63,*) 'm0=',sqrt(m0s),'mtmin=',mtmin1 200 b1=sqrt(p_bm-(bscal_m/en_m)**2) 210 aa=bdotl_m/en_m c 220 write(63,*) 'ellipse plane angle = ',atan(b1/aa)*180./pi 230 kos=aa/sqrt(aa**2+b1**2) 300 lam=w_mass*sqrt((top_m**2-mtmin1**2)/(2.*bdotl_m)) c write(63,*) 'top_m=',top_m,'major _',lam/kos,'minor _',lam circum_m=2*pi*lam/sqrt(kos) cmfrnc_m=circum_m 500 zmt=(top_m**2-m0s)*kos/(2*aa) 1100 al2(1)=b_muo(2)*p_muo(3)-b_muo(3)*p_muo(2) 1101 al2(2)=b_muo(3)*p_muo(1)-b_muo(1)*p_muo(3) 1102 al2(3)=b_muo(1)*p_muo(2)-b_muo(2)*p_muo(1) 1103 al2s=0 1110 do j=1,3 1111 al2s=al2(j)**2+al2s 1112 enddo ! end j c write(63,*)'al2 & al2s',al2(1),al2(2),al2(3),al2s 1120 do j=1,3 1121 al2(j)=al2(j)/sqrt(al2s) 1122 enddo ! end j 1130 al1(1)=(al2(2)*p_muo(3)-al2(3)*p_muo(2))/en_m 1131 al1(2)=(al2(3)*p_muo(1)-al2(1)*p_muo(3))/en_m 1132 al1(3)=(al2(1)*p_muo(2)-al2(2)*p_muo(1))/en_m c write(63,*) 'eta_m','t(x)','t(y)','t(z)' 1200 do k=1,el_grid 1210 alp=(k-1)*2.*pi/(el_grid) eta_m(k)=(k-1)*360./(el_grid) 1220 xp=lam*cos(alp)+(w_mass**2)*(b1/aa)/(2*en_m) 1225 yp=lam*sin(alp) 1230 zp=xp*(b1/aa)+zmt/kos 1300 do j=1,3 1315 top_muo(k,j)=b_muo(j)+(zp+en_m-(w_mass**2/(4.*en_m)))* < p_muo(j)/en_m-xp*al1(j)-yp*al2(j) 1320 enddo ! end j 1390 ts=0 1400 do j=1,3 1410 ts=top_muo(k,j)**2+ts 1420 enddo ! end j c write(63,*) eta_m(k),top_muo(k,1),top_muo(k,2),top_muo(k,3) 2000 enddo ! end k 9009 return end Subroutine cross(x1,x2,fqqb,fglu) c include files c ------------- include 'hadtop_fast.inc' include 'x1x2.inc' include 'ellipse_e.inc' c local declarations c ------------------ real s,shat,pee,kay,that,tlong1,tlong2,e1,e2,theta real et,tep,fglu,fqqb,siggg,sigqq,x1,x2 c Set up kinematics for relative cross sections of q-qbar vs. gluon c fusion. dsig/dthat are given in parton CM c First decide which t is the forward one (with larger positive long- c itudinal momentum if((top_l(3)-top_h(3)).gt.0) then tlong1=top_l(3) e1=top_l(4) tlong2=top_h(3) e2=top_h(4) else tlong1=top_h(3) e1=top_h(4) tlong2=top_l(3) e2=top_l(4) endif c s=(TEV_Energy)**2 shat=x1*x2*s pee=sqrt(shat)/2 if((shat-(2*top_m)**2).lt.0) then write(63,*) 'shat is too small' goto 9000 else endif kay=sqrt((shat/4)-top_m**2) that=top_m**2-x1*sqrt(s)*(e1-tlong1) theta=acos((that-top_m**2+shat/2)/(2*kay*pee)) c Now the two diff cross sections can be written. sigqq=2*(1+cos(theta)**2+4*top_m**2*(sin(theta))**2/shat)/ < (9*shat**2) siggg=(7+9*kay**2*(cos(theta))**2/pee**2)*(1+2*kay**2*top_m**2* < (sin(theta))**2/pee**4-kay**4*(cos(theta))**4/pee**4)/(48* < shat**2*(1-(kay*cos(theta)/pee)**2)**2) c the relative fraction of gluon or qqbar is given by fglu and fqqb fglu=siggg/(siggg+sigqq) fqqb=sigqq/(siggg+sigqq) 9000 return end Subroutine cross_dilep(x1,x2,fqqb,fglu) c include files c ------------- include 'hadtop_fast.inc' include 'x1x2.inc' include 'ellipse_e.inc' include 'ellipse_m.inc' c local declarations c ----------------- real s,shat,pee,kay,that,tlong1,tlong2,e1,e2,theta real et,tep,fglu,fqqb,siggg,sigqq,x1,x2 c Set up kinematics for relative cross sections of q-qbar vs. gluon c fusion. dsig/dthat are given in parton CM c First decide which t is the forward one (with larger positive long- c itudinal momentum if((top_e(3)-top_mu(3)).gt.0) then tlong1=top_e(3) e1=top_e(4) tlong2=top_mu(3) e2=top_mu(4) else tlong1=top_mu(3) e1=top_mu(4) tlong2=top_e(3) e2=top_e(4) endif c s=(TEV_Energy)**2 shat=x1*x2*s pee=sqrt(shat)/2 if((shat-(2*top_m)**2).lt.0) then write(63,*) 'shat is too small' goto 9000 else endif kay=sqrt((shat/4)-top_m**2) that=top_m**2-x1*sqrt(s)*(e1-tlong1) theta=(that-top_m**2+shat/2)/(2*kay*pee) theta=acos((that-top_m**2+shat/2)/(2*kay*pee)) c Now the two diff cross sects can be written. sigqq=2*(1+cos(theta)**2+4*top_m**2*(sin(theta))**2/shat)/ < (9*shat**2) siggg=(7+9*kay**2*(cos(theta))**2/pee**2)*(1+2*kay**2*top_m**2* < (sin(theta))**2/pee**4-kay**4*(cos(theta))**4/pee**4)/(48* < shat**2*(1-(kay*cos(theta)/pee)**2)**2) c the relative fraction of gluon or qqbar is given by fglu and fqqb fglu=siggg/(siggg+sigqq) fqqb=sigqq/(siggg+sigqq) 9000 return end subroutine w_const(jet_iu,jet_id) c calculates set of pairs of jet correction factors to constrain the c 2 jet invariant mass to that of W. c c include files c ------------- include 'hadtop_fast.inc' c local declarations c ------------------ integer iu,id real jet_iu(4),jet_id(4) real jet_u(4),jet_d(4),e_u,e_d,pt_u,pt_d,m_u,m_d,sig_u,sig_d real ptem(4),corr_tem(wm_grid),mtem,pt_tem, prob_u,prob_d integer order,jswitch ! order of x,y;switch 1-2? real fact, mw real phi,pi,psi,c,s,x,y,e1,e2 real coeff1,coeff2,mag1,mag2 real c1,c2,a1,a2,x1,y1,psimin,psimax,x2,y2 data pi /3.14159265359/ c prepare local jet arrays jet_u(1)=jet_iu(1) jet_u(2)=jet_iu(2) jet_u(3)=jet_iu(3) jet_u(4)=jet_iu(4) m_u=sqrt(jet_u(4)**2-jet_u(1)**2-jet_u(2)**2-jet_u(3)**2) c jet_d(1)=jet_id(1) jet_d(2)=jet_id(2) jet_d(3)=jet_id(3) jet_d(4)=jet_id(4) m_d=sqrt(jet_d(4)**2-jet_d(1)**2-jet_d(2)**2-jet_d(3)**2) c c let larger jet mass be jet 2; switch if necessary. if(m_d.lt.m_u) then do jj=1,4 ptem(jj)=jet_u(jj) jet_u(jj)=jet_d(jj) jet_d(jj)=ptem(jj) enddo mtem=m_u m_u=m_d m_d=mtem jswitch=1 else jswitch=0 endif mw=(jet_u(4)+jet_d(4))**2-(jet_u(1)+jet_d(1))**2- < (jet_u(2)+jet_d(2))**2-(jet_u(3)+jet_d(3))**2 mw=sqrt(mw) c write(63,*)'E1= ',jet_u(4),'E2= ',jet_d(4) c write(63,*)'m_u=',m_u,'m_d=',m_d,'mw=',mw c c Now the angle by which an ellipse is rotated in corr_u X corr_d plane c can be between +45o and 0o (not <0 because of switch for m_d0 and phi near 45o c ----------The formula------------------------------- c mw**2=x**2*(c**2*m_u**2+s**2*m_d**2-s*c*fact)+ c < y**2*(s**2*m_u**2+c**2*m_u**2+s*c*fact) c ---------------------------------------------------- coeff1=(c**2*m_u**2+s**2*m_d**2-s*c*fact) coeff2=(s**2*m_u**2+c**2*m_d**2+s*c*fact) c write(63,*)'1st coeff=',(c**2*m_u**2+s**2*m_d**2-s*c*fact) c write(63,*)'2nd coeff=',(s**2*m_u**2+c**2*m_d**2+s*c*fact) c need to decide if((coeff2.lt.0).and.(coeff1.gt.0)) then mag1=w_mass/sqrt(coeff1) mag2=w_mass/sqrt(-coeff2) order=1 else c choose only positive corrections and energy if((coeff1.lt.0).and.(coeff2.gt.0)) then mag2=w_mass/sqrt(coeff2) mag1=w_mass/sqrt(-coeff1) order=2 else print *,'both same sign?' endif endif c write(63,*) 'mag1= ',mag1,' mag2= ',mag2 c get sigmas for the input jet energies e_u=jet_u(4) e_d=jet_d(4) c transverse momenta pt_u=sqrt(jet_u(1)**2+jet_u(2)**2) pt_d=sqrt(jet_d(1)**2+jet_d(2)**2) call msigma(pt_u,sig_u) call msigma(pt_d,sig_d) c write(63,*) ' sig_u= ',sig_u,' sig_d= ',sig_d c scale sig_u and sig_d to full energies sig_u=sig_u*e_u/pt_u sig_d=sig_d*e_d/pt_d c write(63,*) ' sig_eu= ',sig_u,' sig_ed= ',sig_d c find range of psi values for e_u and e_d within chosen range +/-nsigma. c End points are at e_u+nsigma*sig_u or e_d-nsigma*sig_d, whichever is c nearer in to origin, for largest corrected e_u and at e_u-nsigma*sig_u c or ..., for smallest allowed corrected e_u. a1=c**2-(mag2/mag1)**2*s**2 a2=s**2-(mag2/mag1)**2*c**2 c write(63,*) 'a1=',a1,'a2=',a2 do j=1,2 c1=1-(2*j-3)*nsigma*sig_u/e_u c2=1+(2*j-3)*nsigma*sig_d/e_d c write(63,*) 'j=',j,'c1=',c1,'c2=',c2 if(j.eq.1) then if((c1/mag1)**2+a1.lt.0) then x1=0 y1=0 else x1=(c1*c-mag2*s*sqrt((c1/mag1)**2+a1))/a1 y1=mag2*sqrt((x1/mag1)**2+1) endif c write(63,*) 'From c1:c1 out=',x1*c+y1*s,' c2 out=',-x1*s+y1*c if((c2/mag1)**2+a2.lt.0) then x2=0 y2=0 else x2=(-c2*s+mag2*c*sqrt((c2/mag1)**2+a2))/a2 y2=mag2*sqrt((x2/mag1)**2+1) endif c write(63,*) 'From c2:c1 out=',x2*c+y2*s,' c2 out=',-x2*s+y2*c if((abs(c1-x1*c-y1*s).le.c1/1000) < .and.(c2+x1*s-y1*c.lt.0)) then c write(63,*) 'no.1',x1,y1 psimin=alog(x1/mag1+y1/mag2) jmin=1 else if(x2.gt.0) then c write(63,*) 'no.2',x2,y2 psimin=alog(x2/mag1+y2/mag2) jmin=2 else c write(63,*) 'too many zeroes' endif endif else x1=(c1*c-mag2*s*sqrt((c1/mag1)**2+a1))/a1 y1=mag2*sqrt((x1/mag1)**2+1) c write(63,*) 'From c1:c1 out=',x1*c+y1*s,' c2 out=',-x1*s+y1*c x2=(-c2*s+mag2*c*sqrt((c2/mag1)**2+a2))/a2 y2=mag2*sqrt((x2/mag1)**2+1) c write(63,*) 'From c2:c1 out=',x2*c+y2*s,' c2 out=',-x2*s+y2*c if((abs(c1-x1*c-y1*s).le.c1/1000).and. < (c2.gt.-x1*s+y1*c)) then c write(63,*) 'no.3',x1,y1 psimax=alog(x1/mag1+y1/mag2) jmax=1 else c write(63,*) 'no.4',x2,y2 psimax=alog(x2/mag1+y2/mag2) jmax=2 endif endif enddo c write(63,*) 'lower lim on E',jmin,'is psi=',psimin c write(63,*) 'upper lim on E',jmax,'is psi=',psimax c Now I divide the range of allowed psi values into wm_grid points c varying from psimin to psimax. The corrections are the output c along with the probabilities from jet energy gaussians. c Check that the order is y~cosh. If a switch was made, define c corrections oppositely do i=1,wm_grid psi=(i-1)*(psimax-psimin)/(wm_grid-1)+psimin if(order.eq.1) then x=cosh(psi)*mag1 y=sinh(psi)*mag2 else y=cosh(psi)*mag2 x=sinh(psi)*mag1 endif if(jswitch.eq.0) then corr_u(i)=x*c+y*s corr_d(i)=y*c-x*s else corr_d(i)=x*c+y*s corr_u(i)=y*c-x*s endif c write(63,*) 'i= ',i,' corr_u= ',corr_u(i),' corr_d= ',corr_d(i) if((corr_u(i).lt.0).or.(corr_d(i).lt.0)) then write(63,*) 'negative corrections?' goto 101 endif e_u_corr=jet_u(4)*corr_u(i) e_d_corr=jet_d(4)*corr_d(i) c write(63,*)'index= ',i,' e_u_corr= ',e_u_corr, c < ' e_d_corr= ',e_d_corr call gauxx(e_u_corr,e_u,sig_u,prob_u) call gauxx(e_d_corr,e_d,sig_d,prob_d) wm_prob(i)=prob_u*prob_d c write(63,*) prob_u,prob_d, wm_prob(i), jswitch 101 continue enddo c switch the corrections if necessary if(jswitch.eq.1) then do jj=1,wm_grid corr_tem(jj)=corr_u(jj) corr_u(jj)=corr_d(jj) corr_d(jj)=corr_tem(jj) enddo endif return end c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: subroutine had_top(jet_u,jet_d,jet_b) c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c subroutine calculates probabilities for any of the four-momenta c of the top quark decaying into three hadronic jets. Assuming the c measured jet quantities as the central values, each of the three c jets is allowed to vary about its initial values by some epsilon. c The overall probability of a particular (varied) configuration of c jets is a product of Gaussian probabilities for each of the jets c to vary by epsilon from its initial value. c Jets assumed to have originated from a W decaying hadronically into c a pair of quarks have their effective mass constrained to the book c value of MW=80.3 GeV. The third jet, a b candidate, is allowed to c vary within +-NSIG*sigma, where sigma is an error on jet energy c according to the CDF QCD group (subroutine MSIGMA by Naor Wainer). c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c include files c ------------- include 'hadtop_fast.inc' include 'hadpart.inc' c local declarations c ------------------ integer iu,id,ib integer i,j,k c external functions c ------------------ external response_function external degrade c real jet_u(4),jet_d(4),jet_b(4) real jet_u_corr(4) real jet_d_corr(4) real jet_b_corr(4) real pt_u,pt_d,pt_b real top_corr(4) real sig_b, eps_b, prob_b C integer bj_grid,bc_grid,bn_grid C integer wm_grid,wc_grid,wn_grid C integer el_grid,ec_grid,en_grid C C bj_grid=11 C wm_grid=13 C el_grid=37 C bc_grid=(bj_grid+1)/2 C bn_grid=(bj_grid-1)/2 C wc_grid=(wm_grid+1)/2 C wn_grid=(wm_grid-1)/2 C ec_grid=(el_grid+1)/2 C en_grid=(el_grid-1)/2 pt_u=sqrt(jet_u(1)**2+jet_u(2)**2) pt_d=sqrt(jet_d(1)**2+jet_d(2)**2) pt_b=sqrt(jet_b(1)**2+jet_b(2)**2) c c write(63,*) ' jet_u : ',(jet_u(k),k=1,4) c write(63,*) ' jet_d : ',(jet_d(k),k=1,4) c write(63,*) ' jet_b : ',(jet_b(k),k=1,4) c find the energy error for third jet call msigma(pt_b,sig_b) sig_b=sig_b*jet_b(4)/pt_b c write(63,*) pt_b,sig_b c write(63,*) 'm(u,d) uncorrected:',sqrt( c < (jet_u(4)+jet_d(4))**2 c < -(jet_u(1)+jet_d(1))**2 c < -(jet_u(2)+jet_d(2))**2 c < -(jet_u(3)+jet_d(3))**2) c find the corrections for jets originating from W->qq, constrained to MW c call w_const(jet_u,jet_d,u_cor,d_cor) call w_const(jet_u,jet_d) c momenta+energies do i=1,wm_grid ! loop over jets from W->qq decay c write(63,*) 'wm_prob: ',i, wm_prob(i) do j=1,4 jet_u_corr(j)=jet_u(j)*corr_u(i) jet_d_corr(j)=jet_d(j)*corr_d(i) u_jet(i,j)=jet_u_corr(j) d_jet(i,j)=jet_d_corr(j) enddo c test printouts c check that m(u,d) is now at MW c write(63,*) 'm(u,d) corrected:',sqrt( c < (jet_u_corr(4)+jet_d_corr(4))**2 c < -(jet_u_corr(1)+jet_d_corr(1))**2 c < -(jet_u_corr(2)+jet_d_corr(2))**2 c < -(jet_u_corr(3)+jet_d_corr(3))**2), c < ' corr_u=',corr_u(i),' corr_d=',corr_d(i) c build a grid spanning the allowed b-jet energy range do j=1,bj_grid ! b-quark jet eps_b=float(j-bc_grid)/float(bn_grid) eps_b=eps_b*nsigma*sig_b jet_b_corr(4)=jet_b(4)+eps_b c find (Gaussian) probability of such a deviation from the initial b-jet prob_b=response_function(jet_b(4),degrade(jet_b_corr(4))) c call gauxx(jet_b_corr(4),jet_b(4),sig_b,prob_b) eps_b=jet_b_corr(4)/jet_b(4) do k=1,3 jet_b_corr(k)=jet_b(k)*eps_b b_jet(j,k)=jet_b_corr(k) enddo b_jet(j,4)=jet_b_corr(4) c combine the three jets into top four-vector do k=1,3 top_corr(k)=jet_u_corr(k)+jet_d_corr(k)+jet_b_corr(k) top_had(i,j,k)=top_corr(k) enddo top_corr(4)=jet_u_corr(4)+jet_d_corr(4)+jet_b_corr(4) c fill resulting arrays top_had(i,j,4)=top_corr(4) top_mass(i,j)=top_corr(4)**2 < -top_corr(1)**2-top_corr(2)**2-top_corr(3)**2 if(top_mass(i,j).gt.0) top_mass(i,j)=sqrt(top_mass(i,j)) top_prob(i,j)=wm_prob(i)*prob_b c test printouts c write(63,*) (top_corr(k),k=1,4),' mt=',top_mass(i,j) enddo ! j loop b-jet enddo ! i loop MW return end Real Function RESPONSE_FUNCTION(ETRUE,EMEAS) Implicit None Real ROOT2PI, EDEGRADE, SLP1, SLP2, SIG Real XNORM_LOG, PROB_LOG Real ETRUE, EMEAS REAL SIGMA, DEGRADE, SLOPE, FZERO External SIGMA, DEGRADE, DFREQ, SLOPE, FZERO Double Precision DFREQ, CERN, ARG Data ROOT2PI /2.506628/ C EDEGRADE = DEGRADE(ETRUE) C C>>>Find upward (SLP1) and downward (SPL2) slopes: SLP1 = SLOPE(ETRUE,1) SLP2 = SLOPE(ETRUE,2) SIG = SIGMA(ETRUE) C C>>>Evaluate first smeared exponential (upwards): IF (SLP1 .EQ. 0) GOTO 200 XNORM_LOG =(0.5*(SIG/SLP1)**2.)-.5-(EMEAS-EDEGRADE+SLP2/2)/SLP1 ARG = (EDEGRADE - EMEAS - (SLP1 + SLP2)/2)/SIG + (SIG/SLP1) CERN = DFREQ(ARG) If (SLP1 .Gt. 0.0) CERN = 1.0D+00 - CERN If (CERN .Le. 1.0D-37) Then 200 SIG = SQRT(SIG**2 + SLP1**2) C ARG = (EDEGRADE - EMEAS - SLP1 - SLP2)/SIG ARG = (EDEGRADE - EMEAS)/SIG ARG = - (ARG**2.)/2. RESPONSE_FUNCTION = 1.0D+00/(SIG*ROOT2PI) * DEXP(ARG) Else PROB_LOG = XNORM_LOG + DLOG(CERN) RESPONSE_FUNCTION = (1./ABS(SLP1)) * EXP(PROB_LOG) EndIf C C>>>Evaluate second smeared exponential (downwards): IF (SLP2 .EQ. 0) GOTO 201 XNORM_LOG =(0.5*(SIG/SLP2)**2.)-.5-(EMEAS-EDEGRADE+SLP1/2)/SLP2 ARG = (EDEGRADE - EMEAS - (SLP1 + SLP2)/2)/SIG + (SIG/SLP2) CERN = DFREQ(ARG) If (SLP2 .Gt. 0.0) CERN = 1.0D+00 - CERN If (CERN .Le. 1.0D-37) Then 201 SIG = SQRT(SIG**2 + SLP2**2) C ARG = (EDEGRADE - EMEAS - SLP1 - SLP2)/SIG ARG = (EDEGRADE - EMEAS)/SIG ARG = - (ARG**2.)/2. RESPONSE_FUNCTION = (RESPONSE_FUNCTION & + 1.0D+00/(SIG*ROOT2PI) * DEXP(ARG))/2 Else PROB_LOG = XNORM_LOG + DLOG(CERN) RESPONSE_FUNCTION = (RESPONSE_FUNCTION & + (1./ABS(SLP2)) * EXP(PROB_LOG))/2 EndIf RESPONSE_FUNCTION = ABS(1 - FZERO(ETRUE)) * RESPONSE_FUNCTION Return End C C / / / / / / / / / / / / / / / / / / / / C Real Function DEGRADE(ETRUE) Implicit None Real ECM, ETRUE Real A, B, C, D, E, F, a2, OFFSET PARAMETER (A = 0.471) PARAMETER (B = 0.788) PARAMETER (C = 9.831E-04) PARAMETER (D = -3.573E-06) PARAMETER (E = 4.103E-09) PARAMETER (A2 = -0.63) c COMMON/RESPONSE/ECM DATA ECM/1950./ c OFFSET = 0 IF (ABS(ECM-546.) .LT. 10.) OFFSET = A2 If (ETRUE .Ge. 14.5) Then !was 10 8/29/91 DEGRADE = A + B*(ETRUE) + & C*(ETRUE**2) + D*(ETRUE**3) + E*(ETRUE**4) + OFFSET Else DEGRADE = 2.033 + 0.597*ETRUE + 0.00669*(ETRUE**2) + OFFSET EndIf c Return End C C / / / / / / / / / / / / / / / / / / / / C Real Function SLOPE(ETRUE,ISLOPE) COMMON/RESPONSE/ECM C DATA ECM/1950./ C SLOPE = 0 IF (ISLOPE .EQ. 1) THEN SLOPE = 2.164 + 0.01576*ETRUE + (-1.611E-05)*ETRUE**2 IF (ABS(ECM-546.) .LT. 10.) THEN !ECM = 546 SLOPE = SLOPE - 0.48*1.2 !Effect of 546 UE vs. 1950 END IF ELSE IF (ISLOPE .EQ. 2) THEN SLOPE = 0.605 - 0.05046*ETRUE - ( 6.378E-05)*ETRUE**2 IF (SLOPE .GT. -0.2) SLOPE = -0.200 !Low Et has small downward tail END IF RETURN END C C / / / / / / / / / / / / / / / / / / / / C Real Function SIGMA(ETRUE) Implicit None Real ETRUE, EDUM C EDUM = ETRUE If (ETRUE .Lt. 4.) EDUM = 4. SIGMA = 0.320*SQRT(EDUM) + (0.0264)*EDUM + 1.039 - (3.934/EDUM) C Return End C C / / / / / / / / / / / / / / / / / / / / C Real Function FZERO(ETRUE) Implicit None REAL ETRUE INTEGER NEFF DATA NEFF/14/ REAL ET_EFF(14),EFF_JETS(14), EFF_MAX DATA ET_EFF/ & 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, & 12.5, 15.5, 20.5, 26.5/ DATA EFF_JETS/ & 0.079,0.222,0.403,0.548,0.670,0.747,0.796,0.850,0.879,0.895, & 0.934,0.953,0.974,0.981/ LOGICAL LFIRST DATA LFIRST/.TRUE./ REAL XINTERP2 External XINTERP2 C IF (LFIRST) THEN LFIRST = .FALSE. EFF_MAX = XINTERP2(NEFF,ET_EFF,EFF_JETS,ET_EFF(NEFF)) !Max JETS eff END IF FZERO = XINTERP2(NEFF,ET_EFF,EFF_JETS,ETRUE)/EFF_MAX FZERO = 1 - FZERO IF (FZERO .LT. 0) FZERO = 0. C RETURN END C C / / / / / / / / / / / / / / / / / / / / C FUNCTION XINTERP2(NELEMENTS,XAXIS,YAXIS,X) REAL XINTERP2, XAXIS(*),YAXIS(*) IF (X .LT. XAXIS(1)) THEN XINTERP2 = YAXIS(1) RETURN END IF DO IBIN = 1, NELEMENTS-1 IF((X .GE. XAXIS(IBIN)) .AND. (X .LT. XAXIS(IBIN+1))) THEN XINTERP2 = YAXIS(IBIN) & + ((X-XAXIS(IBIN))/(XAXIS(IBIN+1)-XAXIS(IBIN)) ) & * (YAXIS(IBIN+1) - YAXIS(IBIN)) RETURN ENDIF END DO XINTERP2 = YAXIS(NELEMENTS) RETURN END