c csubs.f, version 2.0, November 1, 1995. c This is the workhorse of the program as it returns the wave c function squared. It is too complicated to explain very well, c but should not need to be modified by the user. double precision function corr(kk,mom,kdotr,r) implicit none include 'cparams.f' double precision r,mom,kdotr,eta,pi,arg,krmax,corr0 double precision rmass,wsym1,wsym2,wanti1,wanti2,mass1,mass2 double precision xx,xxprime,xxjj,root2,p1,zk integer kk,iq,npart,l1,l2,jj,i0,lc double precision v,e,e1,e2,gc double precision mc,ddelde double complex cdel1(1000,nka),cdel2(1000,nka), + c1far(nka),c2far(nka),cfar,coul double complex acouly(nka),acfact1(nka),acfact2(nka), + achype(40,nka),achype1(5,nka),achype2(5,nka) double complex ci,cphi1,cphi2,cphis1,cphis2,cphia1,cphia2 integer icheap,ijcheap common /cheapint/icheap,ijcheap double precision cheap double precision jcheap,j1cheap,j2cheap,echeap, + gcheap,rcheap,bcheap,momcheap common /cheapreal/jcheap,j1cheap,j2cheap, + echeap,gcheap,rcheap,bcheap,momcheap common /spwave/cdel1,cdel2,c1far,c2far common /couljunk/acouly,acfact1,acfact2,achype,achype1,achype2 common /corrinfo/rmass,wsym1,wsym2,wanti1,wanti2,iq, + npart,l1,l2,mass1,mass2 pi=3.14159265358979323844d0 cheap=0.d0 ci=(0.d0,1.d0) root2=1.414213562373d0 cphis1=(0.d0,0.d0) cphia1=(0.d0,0.d0) cphis2=(0.d0,0.d0) cphia2=(0.d0,0.d0) if(iq.ne.0) then eta=dfloat(iq)*(rmass/137.036d0)/mom zk=kdotr/mom if(dabs(zk).gt.r) write(66,*) 'zk >r',zk,r cphi1=coul( mom,r,zk,eta,kk) cphi2=coul( mom,r,-zk,eta,kk) else arg=kdotr/197.323d0-2.d0*pi*int(kdotr/(197.323d0*2.d0*pi)) cphi1=cdexp(ci*arg) cphi2=dconjg(cphi1) endif cphis1=(cphi1+cphi2)/root2 cphis2=cphis1 cphia1=(cphi1-cphi2)/root2 cphia2=cphia1 c c Below is a fit to a relativistic Breit Wigner for a p-wave c resonance (lc=1). The correction for putting in the correct c l-value is very minor, except for resonances which overlap the c threshhold. "ddelde" is the derivative of the phase shift c with respect to the c.o.m. energy. The form for the phase shift is c given by: c tan(del) = mom^(2*l+1)*g^2/(e*(mc^2-e^2)) c where g in the above equation is proportional to some coupling c constant, which is dimensionless for the case of two scalars going c through a p-wave resonance. Note that if one uses l=0, that c the correlation function blows up at the origin. This is because c the forumla below assumes that the source is much larger than (1/k) c to use the formula with the corr. function being proportional to the c derivative of the phase shifts. With the p-wave form this is not so c much of a problem. Basically, if the rel. momentum of interest c is small compared to 1/r then one should not use the phase shift c formula, but should use wave function arguments, perhaps based on c eff. range theory. c For a reference on using phase shifts to get correlation functions, c one can read a paper by D.H. Boal, (and I believe Shillcock) which c is somewhere in PRC, but I won't go to the trouble to find it now. c if(icheap.eq.1.and.r/rcheap.lt.10) then e1=dsqrt(mass1**2+mom**2) e2=dsqrt(mass2**2+mom**2) v=(mom/e1)+(mom/e2) e=e1+e2 lc=1 mc=mass1+mass2+echeap cheap=(2*jcheap+1)/((2*j1cheap+1)*(2*j2cheap+1)) cheap=cheap*dexp(-(r**2)/(4*(rcheap**2))) cheap=cheap*(v/(mom*mom)) cheap=cheap*(2*(197.323d0**3))/( 8*dsqrt(pi)*(rcheap**3) ) gc=gcheap*(2*mc*mc/(e*(e+mc)))*(mom/momcheap)**(2*lc+1) ddelde=gc*(e/(e+mc))/( (e-mc)**2+(gc/2.d0)**2 ) ddelde=ddelde*(1-(mc**2-e**2)* + ( 1/(2*e**2) - (2*lc+1)/(2*mom*e*v) )) cheap=cheap*ddelde cheap=cheap*bcheap endif c If there are corrections for strong interactions, add the c change for each partial wave. If npart.eq.0 then there c are no strong int. corrections. if(npart.gt.0) then xx=r*mom/197.323d0 krmax=40.d0*mom/197.323d0 jj=1+int(1000.d0*xx/krmax) if(jj.lt.1001) then p1=kdotr/(r*mom) xxjj=(dfloat(jj)-.5d0)*krmax/1000.d0 if(l1.eq.0) cphis1=cphis1+ + root2*cdel1(jj,kk)/xxjj if(l1.eq.1) cphia1=cphia1+ + root2*p1*3*ci*cdel1(jj,kk)/xxjj if(l1.eq.2) cphis1=cphis1-root2*.5d0*(3*(p1**2)-1.d0) + *5*cdel1(jj,kk)/xxjj c ___________________________ if(npart.gt.1) then if(l2.eq.0) cphis2=cphis2+ + root2*cdel2(jj,kk)/xxjj if(l2.eq.1) cphia2=cphia2+ + root2*p1*3*ci*cdel2(jj,kk)/xxjj if(l2.eq.2) cphis2=cphis2-root2* + .5d0*(3*(p1**2)-1.d0)*5*cdel1(jj,kk)/xxjj endif c ___________________________ else xxprime=xx-2*pi*int(xx/(2.d0*pi)) cfar=cdexp(ci*xxprime-ci*eta*dlog(2.d0*xx)) if(l1.eq.0) cphis1=cphis1-ci*root2*c1far(kk)*cfar/xx if(l1.eq.1) cphia1=cphia1 + -root2*p1*3*c1far(kk)*cfar/xx if(l1.eq.2) cphis1=cphis1+ci*root2*5 + *.5d0*(3*(p1**2)-1.d0)*c1far(kk)*cfar/xx c ___________________________ if(npart.gt.1) then if(l2.eq.0) + cphis2=cphis2-ci*root2*c2far(kk)*cfar/xx if(l2.eq.1) cphia2=cphia2 + -root2*p1*3*c2far(kk)*cfar/xx if(l2.eq.2) cphis2=cphis2+ci*root2*5 + *.5d0*(3*(p1**2)-1.d0)*c1far(kk)*cfar/xx endif c ___________________________ endif endif corr0=wsym1*cphis1*dconjg(cphis1)+wsym2*cphis2*dconjg(cphis2)+ + wanti1*cphia1*dconjg(cphia1)+wanti2*cphia2*dconjg(cphia2) corr=corr0+cheap return end c __________________________________________ subroutine filter1(mass1,px,py,pz,if1) c This should NOT depend on phi! implicit none double precision mass1,px,py,pz,pperp,thetamin,thetamax,theta double precision y,vz,mom(3),mass integer if1 c parameter(pi=3.14159265358979323844d0) c if1=0 c thetamin=(pi/180.d0)*60.d0 c thetamax=(pi/180.d0)*120.d0 c vz=pz/dsqrt(mass1**2+px**2+py**2+pz**2) c y=.5d0*dlog((1.d0+vz)/(1.d0-vz)) c pperp=dsqrt(px**2+py**2) c if(y.gt.-0.5d0.and.y.lt.0.5d0.and.pperp.lt.100.d0) then c if1=1 c else c if1=0 c endif if1=1 return end c ____________________________________________ subroutine filter2(mass1,px1,py1,pz1,e1,mass2, + px2,py2,pz2,e2,ifilter) implicit none double precision mass1,px1,py1,pz1,e1,mass2,px2,py2,pz2,e2, + phimin,phimax,pi,phi,pt,rap,etot integer ifilter,iacc1,iacc2 parameter(pi=3.141592654) c phimin=-pi c phimax=pi c phi1=atan2(py1,px1) c phi2=atan2(py2,px2) c ifilter=0 c if(phi1.gt.phimin.and.phi1.lt.phimax.and.phi2.gt.phimin.and. c + phi2.lt.phimax) ifilter=1 ifilter=1 return end c ____________________________________________ c Here are a couple of sample directional testing routines c If the "ctest" conditions are used, the routines are testing whether c the direction of mom is within a certain angle of the chosen c direction. If "ptest" conditions are used, the test is whether c the orthogonal compontent of the momentum is less than some c certain amount. One can insert whichever condition one wants. c NOTE: If one like to do projections and plot correlation functions c as a function of kpar, kin or kout, then these routine do not c do what you want. These routines merely plot the correlation c functions as a function of dsqrt(kpar**2+kin**2+kout**2) under the c constraint the k is in a certain direction. To plot the c correlation function against kpar,kin, etc., you must make a new c mesh. I don't recommend doing this, but if you must feel free to c call me if you would like assistance. c -------------------------------------------- subroutine accpar(mom,kpar,kin,kout,ia) implicit none double precision mom,kpar,kin,kout,ctest,ptest integer ia ia=0 ctest=dabs(kpar/mom) if(ctest.gt.0.866d0) ia=1 c ptest=dsqrt(kin**2+kout**2) c if(ptest.lt.5) ia=1 return end c ____________________________________________ subroutine accin(mom,kpar,kin,kout,ia) implicit none double precision mom,kpar,kin,kout,ctest,ptest integer ia ia=0 ctest=dabs(kin/mom) if(ctest.gt.0.866d0) ia=1 c ptest=dsqrt(kout**2+kpar**2) c ptest=dabs(kpar) c if(ptest.lt.5) ia=1 return end c ____________________________________________ subroutine accout(mom,kpar,kin,kout,ia) implicit none double precision mom,kpar,kin,kout,ctest,ptest integer ia ia=0 ctest=dabs(kout/mom) c ctest=dabs(kpar/mom) if(ctest.gt.0.866d0) ia=1 c ptest=dsqrt(kpar**2+kin**2) c ptest=dabs(kpar) c if(ptest.lt.5) ia=1 return end c ____________________________________________ c c This routine takes the phase space points and rotates them c and boosts them to the frame necessary to calculate the c wave function. c subroutine transform(m1,px1,py1,pz1,e1,x1,y1,z1,t1, + m2,px2,py2,pz2,e2,x2,y2,z2,t2,mom,r,kdotr) implicit none double precision m1,px1,py1,pz1,e1,x1,y1,z1,t1 double precision m2,px2,py2,pz2,e2,x2,y2,z2,t2,mom,r,kdotr c The variable old is useful in testing the code. double precision k0,kx,ky,kz,k2,x,y,z,t,r2,px,py,pz,p0,p2 double precision pdotr double precision pdotk,pdotp1,pdotp2 integer msame msame=0 if(dabs(m1-m2).lt.0.001d0) msame=1 c px=px1+px2 py=py1+py2 pz=pz1+pz2 p0=e1+e2 p2=p0**2-px**2-py**2-pz**2 c kx=.5d0*(px2-px1) ky=.5d0*(py2-py1) kz=.5d0*(pz2-pz1) k0=.5d0*(e2-e1) if(msame.eq.0) then pdotk=m2**2-m1**2 kx=kx-.5d0*pdotk*px/p2 ky=ky-.5d0*pdotk*py/p2 kz=kz-.5d0*pdotk*pz/p2 k0=k0-.5d0*pdotk*p0/p2 endif mom=dsqrt(dabs(k0**2-kx**2-ky**2-kz**2)) c x=x2-x1 y=y2-y1 z=z2-z1 t=t2-t1 pdotr=p0*t-px*x-py*y-pz*z if(msame.eq.1) then x=x-pdotr*px/p2 y=y-pdotr*py/p2 z=z-pdotr*pz/p2 t=t-pdotr*p0/p2 else pdotp1=p0*e1-px*px1-py*py1-pz*pz1 pdotp2=p0*e2-px*px2-py*py2-pz*pz2 x=x-.5d0*pdotr*(px1/pdotp1+px2/pdotp2) y=y-.5d0*pdotr*(py1/pdotp1+py2/pdotp2) z=z-.5d0*pdotr*(pz1/pdotp1+pz2/pdotp2) t=t-.5d0*pdotr*(e1/pdotp1+e2/pdotp2) endif c r2=t*t-x*x-y*y-z*z kdotr=k0*t-kx*x-ky*y-kz*z c -------- if(r2.gt.0.d0) then write(66,*) 'r2<0',-r2 endif r=dsqrt(dabs(r2)) kdotr=-kdotr if(dabs(kdotr).gt.mom*r) write(66,*) 'kdotr=',kdotr,'mom*r=',mom*r return end c c The following routine finds the 3-components of the momentum c used for directional cuts. c subroutine decompose(m1,px1,py1,pz1,e1, + m2,px2,py2,pz2,e2,kx,ky,kz,mmom,iconv,v_source) implicit none c c iconv=idirconv, the parameter in the main routine which is passed c here. c Choose iconv=1 for convention where kx,ky,kz is defined first in c frame where pz1+pz2=0, then boosted outwards to the frame of the c pair. (in this conv. k^2=Qinv, and the z direction is parallel c to the beam axis) c c Choose iconv=2 for conv. where one boosts to frame of the pair, but c not the outwards boost. (in this conv. the z-direction is parallel c to the beam axis) c c Choose iconv=3 for conv. where one boosts along the z-axis to c a frame specified by vs=constant (in this conv. the z-direction c is along the beam axis) c c Choose iconv=4 for conv. where one boosts along the z-axis to c a frame specified by vs=constant, and chooses the x-axis parallel c to the total momentum (in this conv. the z-axis is not necessarily c parallel to the beam axis, and k^2 does not equal Qinv^2) c c Choose iconv=5 for conv. where one does the same as in convention c 4, only one does one more boost parallel to ptot, such that c k^2=Qinv^2 c integer iconv double precision mmom,m1,px1,py1,pz1,e1,m2,px2,py2,pz2,e2,mom c The variable old is useful in testing the code. double precision old double precision vs,v_source,gamma,v,temp double precision k0,kx,ky,kz,px,py,pz,p0,pperp,ptot,pdotk,p2 integer msame c First boost along the beam axis to the frame of the source c &&&&&& IMPORTANT &&&&&&&&&&&& if(iconv.eq.1.or.iconv.eq.2) vs=(pz1+pz2)/(e1+e2) c If you use a fixed frame set vs below if(iconv.eq.3.or.iconv.eq.4.or.iconv.eq.5) vs=v_source c &&&&&&&&&&&&&&&&&&&&&&&&&&&&&& gamma=1.d0/dsqrt(1-vs**2) msame=0 if(dabs(m1-m2).lt.0.001d0) msame=1 c px=px1+px2 py=py1+py2 pz=pz1+pz2 p0=e1+e2 p2=p0**2-px**2-py**2-pz**2 c kx=.5d0*(px2-px1) ky=.5d0*(py2-py1) kz=.5d0*(pz2-pz1) k0=.5d0*(e2-e1) if(msame.eq.0) then pdotk=m2**2-m1**2 kx=kx-.5d0*pdotk*px/p2 ky=ky-.5d0*pdotk*py/p2 kz=kz-.5d0*pdotk*pz/p2 k0=k0-.5d0*pdotk*p0/p2 endif c pperp=dsqrt(px**2+py**2) c Now boost along the z-axis to the source frame temp=gamma*(k0-vs*kz) kz=gamma*(kz-vs*k0) k0=temp temp=gamma*(p0-vs*pz) pz=gamma*(pz-vs*p0) p0=temp c Now find component out of plane k, (ky=ksidewards,out-of-plane) c and component perp. to z, and in the plane (kx=koutwards) temp=(ky*px-kx*py)/pperp kx=(kx*px+ky*py)/pperp ky=temp if(iconv.eq.1) then c Now boost along axis until in the frame of the pair v=pperp/p0 gamma=1./dsqrt(1.-v**2) kx=gamma*(kx-v*k0) endif c Find k parallel to p, and the in-plane comp. perp. to p if(iconv.eq.4.or.iconv.eq.5) then ptot=dsqrt(pperp**2+pz**2) temp=(kx*pperp+kz*pz)/ptot kz=(kz*pperp-kx*pz)/ptot kx=temp endif if(iconv.eq.5) then v=ptot/p0 gamma=1./dsqrt(1.-v**2) kx=gamma*(kx-v*k0) endif mmom=dsqrt(kx*kx+ky*ky+kz*kz) return end c ______________________________________ subroutine ranrot(px,py,x,y,phi) implicit none double precision px,py,x,y,phi,c,s,px0,x0 c=dcos(phi) s=dsin(phi) px0=px px=c*px0+s*py py=c*py-s*px0 x0=x x=c*x0+s*y y=c*y-s*x0 return end c _______________________________________________ subroutine kkfind(kk,mom,akk,nkmax) implicit none include 'cparams.f' double precision mom,akk(nka),little integer kk,i,nkmax little=-0.00001d0 do 234 i=1,nkmax if(mom.gt.little.and.mom.le.akk(i)) then kk=i go to 235 endif little=akk(i) 234 continue kk=nkmax+1 235 return end c subroutine smear(px,py,pz,psx,psy,psz) implicit none double precision px,py,pz,psx,psy,psz,gauss c psx=px+10.d0*gauss() c psy=py+10.d0*gauss() c psz=pz+10.d0*gauss() psx=px psy=py psz=pz return end c subroutine cheapset() implicit double precision (a-h,o-z) double precision mc,mass1,mass2 integer icheap,ijcheap common /cheapint/icheap,ijcheap double precision jcheap,j1cheap,j2cheap,echeap, + gcheap,rcheap,bcheap,momcheap common /cheapreal/jcheap,j1cheap,j2cheap, + echeap,gcheap,rcheap,bcheap,momcheap if(ijcheap.eq.12) then jcheap=1.5d0 j1cheap=.5d0 j2cheap=0.d0 echeap=1232-938.28-139.58 gcheap=112.d0 bcheap=1.d0 rcheap=1.5d0 mass1=938.28d0 mass2=139.58d0 endif if(ijcheap.eq.11) then jcheap=3.d0 j1cheap=1.d0 j2cheap=0.d0 echeap=.711d0 gcheap=.024d0 rcheap=1.5d0 bcheap=1.d0 mass1=1875.585d0 mass2=3727.323d0 endif if(ijcheap.eq.8) then jcheap=1.d0 j1cheap=0.d0 j2cheap=0.d0 echeap=32.12d0 c Use 4.3 for perfect detector, 7.1 to account for Yufeng's resolution c gcheap=4.43 gcheap=7.1d0 rcheap=1.d0 bcheap=.5d0 mass1=493.646d0 mass2=493.646d0 endif if(ijcheap.eq.13) then gcheap=1.5d0 echeap=1.966d0 rcheap=1.5d0 jcheap=1.5d0 j1cheap=.5d0 j2cheap=0.d0 bcheap=1.d0 mass1=938.28d0 mass2=3726.35d0 endif mc=echeap+mass1+mass2 momcheap=.25d0*mc**2+ + .25d0*((mass1**2-mass2**2)**2)/(mc**2) + -.5d0*(mass1**2+mass2**2) momcheap=dsqrt(momcheap) write(6,*) 'The resonance occurs at a rel. momentum of ',momcheap return end c _______________________________________________ c c This needs to go in the main program c integer*4 IDUM,IDUM2,IY,IV(32) c common/RAN2JUNK/IDUM,IDUM2,IY,IV c double precision ran2 c Set IDUM as a negative integer to initialize c double precision function ran2() implicit none integer IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV real EPS,RNMX,AM parameter (IM1=2147483563,IM2=2147483399,IMM1=(IM1-1)) parameter (IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211) parameter (IR2=3791,NTAB=32,NDIV=(1+IMM1/NTAB)) parameter (EPS=1.2e-7,RNMX=(1.0-EPS),AM=(1.0/IM1)) c integer*4 IDUM integer*4 J integer*4 K integer*4 IDUM2 integer*4 IY integer*4 IV(0:NTAB-1) double precision TEMP common/RAN2JUNK/IDUM,IDUM2,IY,IV c write(6,*) 'IDUM=',IDUM if(IDUM.le.0) then IDUM2=123456789 IY=0 IDUM =-IDUM IDUM2=IDUM do J=NTAB+7,0,-1 K=(IDUM)/IQ1 IDUM=IA1*(IDUM-K*IQ1)-K*IR1 if(IDUM.lt.0) IDUM =IDUM+IM1 if(J.lt.NTAB) IV(J)=IDUM enddo IY=IV(0) c write(6,*) 'Initialization completed, IDUM=', IDUM endif K=(IDUM)/IQ1 IDUM=IA1*(IDUM-K*IQ1)-K*IR1 if(IDUM.lt.0) IDUM=IDUM+IM1 K=IDUM2/IQ2 IDUM2=IA2*(IDUM2-K*IQ2)-K*IR2 if(IDUM2.lt.0) IDUM2=IDUM2+IM2 J=IY/NDIV IY=IV(J)-IDUM2 IV(J)=IDUM if(IY.lt.1) IY=IY+IMM1 TEMP=AM*IY if(TEMP.gt.RNMX) then ran2=RNMX else ran2=TEMP endif return end c double precision function gauss() implicit none double precision x,y,r2,r,c,ran2 100 x=(1.d0-2.d0*ran2()) y=(1.d0-2.d0*ran2()) r2=x**2+y**2 if(r2.gt.1.d0) go to 100 r=dsqrt(r2) c=x/r gauss=c*dsqrt(-2.d0*dlog(r2)) return end