Return to
eekeek.f
CVS log
Up to [HallC] / Poltar
|
Return to
eekeek.f
CVS log
|
Up to [HallC] / Poltar
|
|
File: [HallC] / Poltar / eekeek.f
(download)
Revision: 1.1.1.1 (vendor branch), Wed Oct 22 13:58:52 2003 UTC (20 years, 11 months ago) by jones Branch: poltar, MAIN CVS Tags: start, HEAD Changes since 1.1: +0 -0 lines Import simc poltar |
subroutine eekeek(ss,q22,angl,theta,phi,epsi,sigma_eep)
implicit none
include 'simulate.inc'
real*8 ss,q22,angl,theta,phi,epsi
real*8 w,skc2,skc,q0,q0c,qr,aflx,aflxl,an,x,sx
real*8 ps,qs,as,sstmp,q22tmp
double complex z1a,z2a,z3a,z4a,z7a,z8a,ur,ui
real*8 dsigt00,dsigl00,dsigp00,dsigi00
c real*8 dsigp0,dsigi0,ctt
real*8 zf(2,6),sigma_eep
real*8 qv2,qv,qvc
integer pn,pna(3)
integer i
c real*4 for compatability with CERNLIB routine fint.
real*4 px(3),pa(50)
real*4 fint
c write(21,*) 'eekeek'
sstmp=ss
q22tmp=q22
c write(6,*) ss,q22,angl,theta,phi,epsi,sigma_eep
w=sqrt(ss)*1000.
q22=q22*(-1.d+06)
c write(6,*) w,Mk2,targ.Mrec_struck
skc2=(w**2-Mk2-targ.Mrec_struck**2)**2-4.*Mk2*targ.Mrec_struck**2
c write(6,*) skc2
skc2=max(skc2,0.d0)
skc=sqrt(skc2)/2./w
c write(6,*) -q22,Mp2
q0=-(-q22-w**2+Mp2)/2./Mp
q0c=(-q22+q0*Mp)/w
qr=sqrt(q22)/q0c
qv2=q22+q0**2
qv=sqrt(qv2)
qvc=Mp/w*qv
C ctt not used, and causes floating exception if skc=0 (i.e. skc2<0)
C ctt=pi/skc/qvc*1.e+06
aflx=skc/2./w/(w**2-Mp2)*(hbarc)**2*10000.
aflxl=aflx*qr**2
an=angl*180./pi
c write(21,*) q0,q0c,qr**2*epsi
x=cos(angl)
sx=sin(angl)
pn=3
px(1)=ss
px(2)=q22/1.e+06
px(3)=an
pna(1)=10
pna(2)=11
pna(3)=19
ps=2.6
qs=0.0
as=0.0
do i=1,10
pa(i)=ps
ps=ps+0.3
enddo
do i=11,21
pa(i)=qs
qs=qs+0.2
enddo
do i=22,40
pa(i)=as
as=as+10.
enddo
c write(6,*) ps,qs,as
c write(6,*) pn
c write(21,*) px(1),px(2),px(3)
c write(6,*) pna(1),pna(2),pna(3)
c do i=1,40
c write(6,*) pa(i)
c enddo
zf(1,1)=dble(fint(pn,px,pna,pa,zrff1))
c write(6,*) zf(1,1)
zf(1,2)=dble(fint(pn,px,pna,pa,zrff2))
zf(1,3)=dble(fint(pn,px,pna,pa,zrff3))
zf(1,4)=dble(fint(pn,px,pna,pa,zrff4))
zf(1,5)=dble(fint(pn,px,pna,pa,zrff5))
zf(1,6)=dble(fint(pn,px,pna,pa,zrff6))
zf(2,1)=dble(fint(pn,px,pna,pa,ziff1))
zf(2,2)=dble(fint(pn,px,pna,pa,ziff2))
zf(2,3)=dble(fint(pn,px,pna,pa,ziff3))
zf(2,4)=dble(fint(pn,px,pna,pa,ziff4))
zf(2,5)=dble(fint(pn,px,pna,pa,ziff5))
zf(2,6)=dble(fint(pn,px,pna,pa,ziff6))
c do i=1,6
c write(21,*) zf(1,i),zf(2,i)
c enddo
ur=(1.D0,0.D0)
ui=(0.D0,1.D0)
z1a=zf(1,1)*ur+zf(2,1)*ui
z2a=zf(1,2)*ur+zf(2,2)*ui
z3a=zf(1,3)*ur+zf(2,3)*ui
z4a=zf(1,4)*ur+zf(2,4)*ui
z7a=zf(1,5)*ur+zf(2,5)*ui
z8a=zf(1,6)*ur+zf(2,6)*ui
c write(21,*) z1a,z2a,z3a,z4a
c write(21,*) z7a,z8a
c
c calculate free cross section using Saghai's form.
c
dsigt00=aflx*(abs(z1a)**2+abs(z2a)**2
& +2.*real(conjg(z1a)*z2a)*x+0.5*sx**2*(abs(z3a)**2
& +abs(z4a)**2+2.*real(conjg(z1a)*z4a-conjg(z2a)*z3a
& +conjg(z3a)*z4a*x)))
dsigl00=aflxl*epsi*(abs(z7a)**2+abs(z8a)**2
& +2.*real(conjg(z7a)*z8a)*x)
dsigp00=aflx*epsi*(sin(theta))**2*cos(2.*phi)*
& (0.5*abs(z3a)**2+0.5*abs(z4a)**2+real(conjg(z1a)*z4a-
& conjg(z2a)*z3a+conjg(z3a)*z4a*x))
dsigi00=aflx*sqrt(2.*qr**2*epsi*(1.+epsi))*
& sin(theta)*cos(phi)*real((z7a)*(conjg(z3a)-conjg(z2a)+
& conjg(z4a)*x)+z8a*(conjg(z1a)+conjg(z3a)*x+conjg(z4a)))
C dsigp0=aflx*epsi*(sin(theta))**2*
C & (0.5*abs(z3a)**2+0.5*abs(z4a)**2+real(conjg(z1a)*z4a-
C & conjg(z2a)*z3a+conjg(z3a)*z4a*x))
C dsigi0=aflx*sqrt(2.*qr**2*epsi*(1.+epsi))*
C & sin(theta)*real(z7a*(conjg(z3a)-conjg(z2a)+
C & conjg(z4a)*x)+z8a*(conjg(z1a)+conjg(z3a)*x+conjg(z4a)))
sigma_eep=dsigt00+dsigl00+dsigp00+dsigi00
c write(21,*) 'dsig',dsigt00,dsigl00/epsi,sigma_eep
c write(21,*) 'dsig',dsigt00*ctt,dsigl00/epsi*ctt,sigma_eep*ctt
c write(21,*) dsigp00,dsigi00,sigma_eep
c write(21,*) dsigp0,dsigi0,sigma_eep
ss=sstmp
q22=q22tmp
return
end
subroutine eekeeks(ss,q22,angl,theta,phi,epsi,sigma_eep)
include 'simulate.inc'
real*8 ss,q22,angl,theta,phi,epsi
real*8 w,skc2,skc,q0,q0c,qr,aflx,aflxl,an,x,sx
real*8 as,sstmp,q22tmp
double complex z1a,z2a,z3a,z4a,z7a,z8a,ur,ui
real*8 dsigt00,dsigl00,dsigp00,dsigi00
real*8 zf(2,6),sigma_eep
real*8 qv2,qv,qvc,ctt
integer pn,pna(3),i
c real*4 for compatability with CERNLIB routine fint.
real*4 px(3),pa(50)
real*4 fint
sstmp=ss
q22tmp=q22
w=sqrt(ss)*1000.
q22=q22*(-1.d+06)
skc2=(w**2-Mk2-targ.Mrec_struck**2)**2-4.*Mk2*targ.Mrec_struck**2
skc2=max(skc2,0.d0)
skc=sqrt(skc2)/2./w
q0=-(-q22-w**2+Mp2)/2./Mp
q0c=(-q22+q0*Mp)/w
qr=sqrt(q22)/q0c
qv2=q22+q0**2
qv=sqrt(qv2)
qvc=Mp/w*qv
C ctt not used, and causes floating exception if skc=0 (i.e. skc2<0)
C ctt=pi/skc/qvc*1.e+06
aflx=skc/2./w/(w**2-Mp2)*(hbarc)**2*10000.
aflxl=aflx*qr**2
an=angl*180./pi
x=cos(angl)
sx=sin(angl)
pn=3
px(1)=ss
px(2)=q22/1.e+06
px(3)=an
c write(12,*) px(1),px(2),px(3)
pna(1)=20
pna(2)=10
pna(3)=19
as=0.0
pa(1)=2.851
pa(2)=2.898
pa(3)=2.945
pa(4)=2.991
pa(5)=3.038
pa(6)=3.085
pa(7)=3.132
pa(8)=3.320
pa(9)=3.507
pa(10)=3.695
pa(11)=3.883
pa(12)=4.070
pa(13)=4.258
pa(14)=4.446
pa(15)=4.633
pa(16)=4.821
pa(17)=5.009
pa(18)=5.196
pa(19)=5.384
pa(20)=5.572
pa(21)=0.0
pa(22)=0.250
pa(23)=0.376
pa(24)=0.520
pa(25)=0.750
pa(26)=1.000
pa(27)=1.250
pa(28)=1.500
pa(29)=1.750
pa(30)=2.000
do i=31,49
pa(i)=as
as=as+10.
enddo
zf(1,1)=dble(fint(pn,px,pna,pa,zsrff1))
zf(1,2)=dble(fint(pn,px,pna,pa,zsrff2))
zf(1,3)=dble(fint(pn,px,pna,pa,zsrff3))
zf(1,4)=dble(fint(pn,px,pna,pa,zsrff4))
zf(1,5)=dble(fint(pn,px,pna,pa,zsrff5))
zf(1,6)=dble(fint(pn,px,pna,pa,zsrff6))
zf(2,1)=dble(fint(pn,px,pna,pa,zsiff1))
zf(2,2)=dble(fint(pn,px,pna,pa,zsiff2))
zf(2,3)=dble(fint(pn,px,pna,pa,zsiff3))
zf(2,4)=dble(fint(pn,px,pna,pa,zsiff4))
zf(2,5)=dble(fint(pn,px,pna,pa,zsiff5))
zf(2,6)=dble(fint(pn,px,pna,pa,zsiff6))
ur=(1.D0,0.D0)
ui=(0.D0,1.D0)
z1a=zf(1,1)*ur+zf(2,1)*ui
z2a=zf(1,2)*ur+zf(2,2)*ui
z3a=zf(1,3)*ur+zf(2,3)*ui
z4a=zf(1,4)*ur+zf(2,4)*ui
z7a=zf(1,5)*ur+zf(2,5)*ui
z8a=zf(1,6)*ur+zf(2,6)*ui
c
c calculate free cross section using Saghai's form.
c
dsigt00=aflx*(abs(z1a)**2+abs(z2a)**2
& +2.*real(conjg(z1a)*z2a)*x+0.5*sx**2*(abs(z3a)**2
& +abs(z4a)**2+2.*real(conjg(z1a)*z4a-conjg(z2a)*z3a
& +conjg(z3a)*z4a*x)))
dsigl00=aflxl*epsi*(abs(z7a)**2+abs(z8a)**2
& +2.*real(conjg(z7a)*z8a)*x)
dsigp00=aflx*epsi*(sin(theta))**2*cos(2.*phi)*
& (0.5*abs(z3a)**2+0.5*abs(z4a)**2+real(conjg(z1a)*z4a-
& conjg(z2a)*z3a+conjg(z3a)*z4a*x))
dsigi00=aflx*sqrt(2.*qr**2*epsi*(1.+epsi))*
& sin(theta)*cos(phi)*real((z7a)*(conjg(z3a)-conjg(z2a)+
& conjg(z4a)*x)+z8a*(conjg(z1a)+conjg(z3a)*x+conjg(z4a)))
sigma_eep=dsigt00+dsigl00+dsigp00+dsigi00
c write(12,*) sigma_eep,dsigt00,dsigl00
c write(12,*) dsigp00,dsigi00
ss=sstmp
q22=q22tmp
return
end
subroutine phspwght(delta,theta,phi,weight)
include 'simulate.inc'
real*8 delta,theta,phi
real*8 weight
integer pn,pna(3),i
real*8 ps,ts,ds
c real*4 for compatability with CERNLIB routine fint.
real*4 px(3),pa(78)
real*4 fint
if(electron_arm.eq.1 .and. hadron_arm.eq.2)then !e- in HMS, K in SOS
pn=2
px(1)=theta
px(2)=phi
pna(1)=20
pna(2)=50
ts=-0.03325
ps=-0.0735
do i=1,20
pa(i)=ts
ts=ts+0.0035
enddo
do i=21,70
pa(i)=ps
ps=ps+0.003
enddo
weight=dble(fint(pn,px,pna,pa,weightc))
else if (electron_arm.eq.2 .and. hadron_arm.eq.1) then !e- in SOS,K in HmS
pn=3
px(1)=delta
px(2)=theta
px(3)=phi
pna(1)=8
pna(2)=40
pna(3)=30
ts=-0.063375
ps=-0.0435
ds=-17.5
do i=1,8
pa(i)=ds
ds=ds+5.
enddo
do i=9,48
pa(i)=ts
ts=ts+0.00325
enddo
do i=49,78
pa(i)=ps
ps=ps+0.003
enddo
weight=dble(fint(pn,px,pna,pa,weightd))
else
write(6,*) 'electron_arm=',electron_arm,' and hadron_arm=',hadron_arm
write(6,*) 'eekeek.f has a phase space factor that is only defined for'
write(6,*) 'hms&sos case. Need to update for other spectrometers.'
stop
endif
weight=max(weight,0.01D00)
weight=max(100.D00/weight,1.0D00)
return
end
| Analyzer/Replay: Mark Jones, Documents: Stephen Wood |
Powered by ViewCVS 0.9.2-cvsgraph-1.4.0 |