1 gaskelld 1.1 	real*8 function peeK(vertex,main,survivalprob)
  2              
  3              *
  4              * Purpose:
  5              * This routine calculates N(e,e'K)Y cross sections.
  6              *
  7              * variables:
  8              *
  9              *   output:
 10              *	peeK	!d5sigma/(dE_e'*dOmega_e'*Omega_K) - fm^2/MeV/sr^2 ???
 11              *
 12              	implicit none
 13              	include 'simulate.inc'
 14              
 15              	type(event_main):: main
 16              	type(event):: vertex
 17              
 18              	real*8 S,phi
 19              	real*8 sigma_eek,sigma_saghai
 20              	real*8 k_eq,gtpr,fac
 21              	real*8 tfcos,tfsin
 22 gaskelld 1.1 	real*8 pathlen,zaero,betak,gammak,p_kaon
 23              	real*8 survivalprob
 24              
 25              ! Variables calculated in transformation to gamma-NUCLEON center of mass.
 26              	real*8 gstar,bstar,bstarx,bstary,bstarz		!beta of boost to C.M.
 27              	real*8 nustar,qstar,qstarx,qstary,qstarz	!q in C.M.
 28              	real*8 ekcm,pkcm,pkcmx,pkcmy,pkcmz		!p_hadron in C.M.
 29              	real*8 ebeamcm,pbeamcm,pbeamcmx,pbeamcmy,pbeamcmz !p_beam in C.M.
 30              	real*8 etarcm,ptarcm,ptarcmx,ptarcmy,ptarcmz	!p_fermi in C.M.
 31              	real*8 thetacm,phicm,phiqn,jacobian,jac_old
 32              
 33              	real*8 sig_factorized
 34              
 35              * Calculate velocity of PHOTON-NUCLEON C.M. system in the lab frame. Use beta
 36              * and gamma of the cm system (bstar and gstar) to transform particles into
 37              * c.m. frame.  Define z along the direction of q, and x to be along the
 38              * direction of the pion momentum perpendicular to q.
 39              
 40              	call transform_to_cm(vertex,main,
 41                   &		gstar,bstar,bstarx,bstary,bstarz,
 42                   &		nustar,qstar,qstarx,qstary,qstarz,
 43 gaskelld 1.1      &		ekcm,pkcm,pkcmx,pkcmy,pkcmz,
 44                   &		ebeamcm,pbeamcm,pbeamcmx,pbeamcmy,pbeamcmz,
 45                   &		etarcm,ptarcm,ptarcmx,ptarcmy,ptarcmz,
 46                   &		thetacm,phicm,phiqn,jacobian,jac_old)
 47              
 48              * Jacobian goes from dt dphi_cm --> dOmega_lab
 49              * Cross section is in terms of dOmega_cm (not dt dphi_cm), so need
 50              * to include dOmega_cm --> dt dphi_cm jacobian (which is 1/2.*pk_cm*q_cm)
 51              	jacobian = jacobian/(2.*pkcm*qstar)
 52              	jac_old = jac_old/(2.*pkcm*qstar)
 53              
 54              	main%thetacm = thetacm
 55              	main%phicm = phicm
 56              	main%pcm = pkcm
 57              	main%davejac = jacobian
 58              	main%johnjac = jac_old          !approx. assuming collinear boost.
 59              !	write (6,*) jacobian,jac_old,100.*(jacobian-jac_old)/jacobian,'%'
 60              
 61              
 62              * calculate some kinematical variables
 63              * 'f' and 'fer' indicate fermi momenta. 'star' or 'cm' indicate CM system
 64 gaskelld 1.1 
 65              	tfcos = pferx*vertex%uq%x+pfery*vertex%uq%y+pferz*vertex%uq%z
 66              	if(tfcos-1..gt.0..and.tfcos-1..lt.1.d-8)tfcos=1.0 
 67              	tfsin = sqrt(1.-tfcos**2)
 68              
 69              	s = (vertex%nu+efer)**2-(vertex%q+pfer*tfcos)**2-(pfer*tfsin)**2
 70              	main%wcm = sqrt(s) 
 71              
 72              
 73              *******************************************************************************
 74              * Get cross section in photon-nucleon center of mass.  sigcm1 is Saghai
 75              * model (in eekeek.f). sigcm2 is factorized model (Doug Koltenuk).
 76              
 77              
 78              *NOTE: wants odd sign/units for s,q2,etc...  Also, theta is calculated
 79              * above, but phi is not, SO WE ARE ALWAYS CALCULATING FOR PHI=0!!
 80              
 81              	if (targ%Mrec_struck.lt.1150.) then	!lambda production.
 82              	  call eekeek(s/1.d6,vertex%q2,main%thetacm,main%theta_pq,phi,main%epsilon,sigma_saghai)
 83              	else
 84              	  call eekeeks(s/1.d6,vertex%q2,main%thetacm,main%theta_pq,phi,main%epsilon,sigma_saghai)
 85 gaskelld 1.1 	endif
 86              	ntup%sigcm1 = sigma_saghai
 87              
 88              
 89              * Factorization model (CURRENTLY SET UP FOR HYDROGEN ONLY!!!)
 90              
 91              	ntup%sigcm2 = sig_factorized(vertex%q2,main%wcm,main%t,pkcm,targ%Mrec_struck)
 92              
 93              * Choose the cross section model to use by default.
 94              !	sigma_eek = ntup.sigcm1		!Saghai
 95              	sigma_eek = ntup%sigcm2		!Factorized model
 96              
 97              	ntup%sigcm = sigma_eek
 98              
 99              
100              * Virtual photon to electron beam flux conversion factor
101              * DJG,2000: Replace targ.Mtar_struck in denominator of gammaflux with more
102              * general efer-pfer*tfcos, for pfer =0 this reverts to old form
103              !       k_eq = (s-targ%Mtar_struck**2)/2./(efer-pfer*tfcos)
104              
105              * JRA,2001: Go back to original version - more consistent with phase space used
106 gaskelld 1.1 * in the subroutine (according to DJG - see gaskell_model.ps)
107                      k_eq = (main%wcm**2-targ%Mtar_struck**2)/2./targ%Mtar_struck
108              
109              
110              * 'sig' is two-fold C.M. cross section: d2sigma/Omega_cm [ub/sr].
111              * Convert from dt dOmega_cm --> dOmega_lab using 'jacobian' [ub/sr]
112              * Convert to 5-fold by multiplying by flux factor, gtpr [1/MeV]
113              * to give d5sigma/dOmega_pi/dOmega_e/dE_e [ub/Mev/sr]
114              * Note that 'jacobian' above is the full jacobian from generate_cm,
115              * (dt dphi_cm --> dOmega_lab) times 2*pkcm*qcm (dOmega_cm --> dt dphi_cm).
116              *
117              * Note that there is an additional factor 'fac' included with gtpr.   This
118              * takes into account pieces in the flux factor that are neglected (=1) in
119              * colinear collisions.  The flux factor is |v_1-v_2| * 2E_1 * 2E_2.
120              * For a stationary target, v_2=0 and so velocity term is v_1=1 (electron
121              * beam), and E_2=M_2.  For collinear boost, the flux factor can be expressed
122              * in a way that is lorenz invariant, and so can be used for lab or C.M.
123              * For a NON-COLLINEAR boost, there are two changes.  First, the |v| term
124              * becomes 1 - (z component of pfer)/efer.  Second, E_2 isn't just the mass,
125              * it becomes E_fermi, so we have to remove targ.Mtar_struck (which is used
126              * for E_2 by default) and replace it with efer.  Since the flux factor
127 gaskelld 1.1 * comes in the denominator, we replace the usual flux factor (gtpr) with
128              * gtpr*fac, where fac = 1/ ( (1-pfer_z/efer)* (efer/mtar_struck) ).
129              
130              
131                      fac = 1./(1.-pferz*pfer/efer) * targ%Mtar_struck/efer
132              	gtpr = alpha/2./(pi**2)*vertex%e%E/vertex%Ein*k_eq/vertex%q2/(1.-main%epsilon)
133              
134              	peeK = sigma_eek*jacobian*(gtpr*fac) !ub/MeV^2/rad-->ub/sr-->ub/MeV/sr
135              
136              
137              c If doing_decay=.false., generate survival probability for main.weight.
138              c main.FP.p.path is dist. to back of detector, want decay prob. at aerogel.
139              C NOTE THAT ZAERO IS TAKEN WITH RESPECT TO THE POSITION AT WHICH PATHLEN
140              C IS CALCULATED (USUALLY THE BACK OF THE CALORIMETER).  IF THE DRIFTS IN
141              C MC_SOS_HUT ARE CHANGED, THEN THE STARTING POINT MAY BE DIFFERENT AND THESE
142              C NUMBERS MAY BE WRONG.  AERO BETWEEN S2Y AND S2X IN SOS.
143              
144              C Beta/Gamma for decay need to use momentum after radiation/eloss, not vertex
145              C momentum.  Get particle momentum from main%SP%p%delta
146              
147              	if (.not.doing_decay) then
148 gaskelld 1.1 	  if (hadron_arm.eq.1) then
149              	    zaero = 0.			!no aerogel yet, use full length.
150              	  else if (hadron_arm.eq.2) then
151              *	    zaero = -76.		!aero at 270cm,last project=346(cal).
152              	    zaero = -82.8		!From Rick: aero at 263.2cm,last project=346(cal).
153              	  else if (hadron_arm.eq.3) then
154              	    zaero = -183.		!aero at 130cm,last project=313(S2)
155              	  else if (hadron_arm.eq.4) then
156              	    zaero = -183.
157              	  endif
158              	  pathlen = main%FP%p%path + zaero*(1+main%FP%p%dx**2+main%FP%p%dy**2)
159              	  p_kaon = spec%p%P*(1.+main%SP%p%delta/100.)
160              	  betak = spec%p%P/sqrt(spec%p%P**2+Mh2)
161              	  gammak = 1./sqrt(1.-betak**2)
162              	  survivalprob = 1./exp(pathlen/(ctau*betak*gammak))
163              	  decdist = survivalprob		!decdist in ntuple
164              	endif
165              
166              	return
167              	end
168              
169 gaskelld 1.1 
170              
171              
172              	real*8 function sig_factorized(q2,w,t,pk,mrec)
173              
174              * Purpose:
175              * This routine calculates p(e,e'K+)Lambda cross sections using a
176              * factorized cross section model from Doug Koltenuk's thesis.
177              *
178              * Cross section (at theta_cm=0, i.e. t=tmin) is F(W)*F(Q^2).
179              * The t-dependance at fixed Q^2,W is F(t)~exp(-A*t),
180              * so the cross section is F(W)*F(Q^2)*(F(t)/F(tmin))=F(W)*F(Q^2)*F(t-tmin)
181              *
182              * The model requires W, Q^2, t, and pkcm.
183              
184              	real*8 q2,w,t,pk,mrec			!all in Mev,MeV**2
185              	real*8 nu,q,nucm,qcm,tmin,pktest
186              	real*8 w2val,q2val,tval,pkval,tminval	!Vars. used in model (GeV)
187              	real*8 fact_w,fact_q,fact_t
188              
189              	include 'constants.inc'
190 gaskelld 1.1 
191              ! Initialize some stuff.  Start with intermediate variables, all in MeV.
192              
193              	nu = ( w**2 + q2 - mp**2 )/2./Mp
194              	q = sqrt(q2+nu**2)
195              	qcm = q*(Mp/W)
196              	nucm = sqrt(qcm**2-q2)
197              	tmin = -1.*(Mk2 - q2 - 2*nucm*sqrt(pk**2+mk2) + 2*qcm*pk )
198              
199              ! Check center of mass pk, since we can get it from w2
200              	pktest = sqrt ( (w**2+mk2-mrec**2)**2/4./w**2 - mk2 )
201              	if (abs((pktest-pk)/pk).gt.0.001) then
202              	  write(6,*) 'Kaon C.M. momentum passed to sig_factorized does not agree with'
203              	  write(6,*) 'the value calculated from W'
204              	  write(6,*) 'Passed,calculated=',pk,pktest
205              	endif
206              
207              ! Parameters used by the model, all converted to GeV.
208              	q2val = q2/1.d6
209              	w2val = w**2/1.d6
210              	pkval = pk/1000.
211 gaskelld 1.1 	tval = t/1.d6
212              	tminval = tmin/1.d6
213              
214              	if (mrec.lt.1150.) then	!lambda production
215              	  fact_q = 1./(q2val+2.67)**2
216              	  fact_t = exp(-2.1*(tval-tminval))
217              	  if (w2val.ne.0) then			!=0 for central event!!!
218              	   fact_w = 0.959*4.1959*pkval/(sqrt(w2val)*(w2val-0.93827**2))
219              	   fact_w = fact_w + (0.18*1.72**2*0.10**2) /
220                   >		( (w2val-1.72**2)**2 + 1.72**2*0.10**2 )
221              	  endif
222              	else
223              	  fact_q = 1./(q2val+0.79)**2
224              	  fact_t = exp(-1.0*(tval-tminval))
225              	  if (w2val.ne.0) then			!=0 for central event!!!
226              	   fact_w = 0.959*4.1959*pkval/(sqrt(w2val)*(w2val-0.93827**2))
227              
228              !	   fact_w = fact_w + (0.18*1.72**2*0.10**2) /
229              !     >		( (w2val-1.72**2)**2 + 1.72**2*0.10**2 )
230              	  endif
231              	endif
232 gaskelld 1.1 
233              	sig_factorized = fact_q*fact_t*fact_w
234              
235              	return
236              	end
237              
238              
239              
240              	subroutine eekeek(ss,q22,angl,theta,phi,epsi,sigma_eep)
241              
242              	implicit none
243              
244              	include 'simulate.inc'
245              
246              	real*8 ss,q22,angl,theta,phi,epsi
247              	real*8 w,skc2,skc,q0,q0c,qr,aflx,aflxl,an,x,sx
248              	real*8 ps,qs,as
249              	double complex z1a,z2a,z3a,z4a,z7a,z8a,ur,ui
250              	real*8 dsigt00,dsigl00,dsigp00,dsigi00
251              c	real*8 dsigp0,dsigi0,ctt
252              	real*8 zf(2,6),sigma_eep
253 gaskelld 1.1 	real*8 qv2,qv,qvc
254              	integer pn,pna(3)
255              	integer i
256              
257              c real*4 for compatability with CERNLIB routine fint.
258              	real*4 px(3),pa(50)
259              
260              	real*4 fint
261              
262              	w=sqrt(ss)*1000.
263              	skc2=(w**2-Mk2-targ%Mrec_struck**2)**2-4.*Mk2*targ%Mrec_struck**2
264              	skc2=max(skc2,0.d0)
265              	skc=sqrt(skc2)/2./w
266              	q0=-(-q22-w**2+Mp2)/2./Mp
267              	q0c=(-q22+q0*Mp)/w
268              	qr=sqrt(q22)/q0c
269              	qv2=q22+q0**2
270              	qv=sqrt(qv2)
271              	qvc=Mp/w*qv
272              C ctt not used, and causes floating exception if skc=0 (i.e. skc2<0)
273              C	ctt=pi/skc/qvc*1.d+06
274 gaskelld 1.1 	aflx=skc/2./w/(w**2-Mp2)*(hbarc)**2*10000.
275              	aflxl=aflx*qr**2
276              	an=angl*180./pi
277              	x=cos(angl)
278              	sx=sin(angl)
279              	pn=3
280              	px(1)=ss
281              	px(2)=q22/1.d+06
282              	px(3)=an
283              	pna(1)=10
284              	pna(2)=11
285              	pna(3)=19
286              	ps=2.6
287              	qs=0.0
288              	as=0.0
289              	do i=1,10
290              	  pa(i)=ps
291              	  ps=ps+0.3
292              	enddo
293              	do i=11,21
294              	  pa(i)=qs
295 gaskelld 1.1 	  qs=qs+0.2
296              	enddo
297              	do i=22,40
298              	  pa(i)=as
299              	  as=as+10.
300              	enddo
301              
302              	zf(1,1)=dble(fint(pn,px,pna,pa,zrff1))
303              	zf(1,2)=dble(fint(pn,px,pna,pa,zrff2))
304              	zf(1,3)=dble(fint(pn,px,pna,pa,zrff3))
305              	zf(1,4)=dble(fint(pn,px,pna,pa,zrff4))
306              	zf(1,5)=dble(fint(pn,px,pna,pa,zrff5))
307              	zf(1,6)=dble(fint(pn,px,pna,pa,zrff6))
308              	zf(2,1)=dble(fint(pn,px,pna,pa,ziff1))
309              	zf(2,2)=dble(fint(pn,px,pna,pa,ziff2))
310              	zf(2,3)=dble(fint(pn,px,pna,pa,ziff3))
311              	zf(2,4)=dble(fint(pn,px,pna,pa,ziff4))
312              	zf(2,5)=dble(fint(pn,px,pna,pa,ziff5))
313              	zf(2,6)=dble(fint(pn,px,pna,pa,ziff6))
314              
315              	ur=(1.D0,0.D0)
316 gaskelld 1.1 	ui=(0.D0,1.D0)  
317              
318              	z1a=zf(1,1)*ur+zf(2,1)*ui
319              	z2a=zf(1,2)*ur+zf(2,2)*ui
320              	z3a=zf(1,3)*ur+zf(2,3)*ui
321              	z4a=zf(1,4)*ur+zf(2,4)*ui
322              	z7a=zf(1,5)*ur+zf(2,5)*ui
323              	z8a=zf(1,6)*ur+zf(2,6)*ui
324              c
325              c	calculate free cross section using Saghai's form.
326              c
327              	dsigt00=aflx*(abs(z1a)**2+abs(z2a)**2
328                   &		+2.*real(conjg(z1a)*z2a)*x+0.5*sx**2*(abs(z3a)**2
329                   &		+abs(z4a)**2+2.*real(conjg(z1a)*z4a-conjg(z2a)*z3a
330                   &		+conjg(z3a)*z4a*x)))
331              	dsigl00=aflxl*epsi*(abs(z7a)**2+abs(z8a)**2
332                   &		+2.*real(conjg(z7a)*z8a)*x)
333              	dsigp00=aflx*epsi*(sin(theta))**2*cos(2.*phi)*
334                   &		(0.5*abs(z3a)**2+0.5*abs(z4a)**2+real(conjg(z1a)*z4a-
335                   &		conjg(z2a)*z3a+conjg(z3a)*z4a*x))
336              	dsigi00=aflx*sqrt(2.*qr**2*epsi*(1.+epsi))*
337 gaskelld 1.1      &		sin(theta)*cos(phi)*real((z7a)*(conjg(z3a)-conjg(z2a)+
338                   &		conjg(z4a)*x)+z8a*(conjg(z1a)+conjg(z3a)*x+conjg(z4a)))
339              C	dsigp0=aflx*epsi*(sin(theta))**2*
340              C     &		(0.5*abs(z3a)**2+0.5*abs(z4a)**2+real(conjg(z1a)*z4a-
341              C     &		conjg(z2a)*z3a+conjg(z3a)*z4a*x))
342              C	dsigi0=aflx*sqrt(2.*qr**2*epsi*(1.+epsi))*
343              C     &		sin(theta)*real(z7a*(conjg(z3a)-conjg(z2a)+
344              C     &		conjg(z4a)*x)+z8a*(conjg(z1a)+conjg(z3a)*x+conjg(z4a)))
345              
346              	sigma_eep=dsigt00+dsigl00+dsigp00+dsigi00
347              c	write(21,*) 'dsig',dsigt00,dsigl00/epsi,sigma_eep
348              c	write(21,*) 'dsig',dsigt00*ctt,dsigl00/epsi*ctt,sigma_eep*ctt
349              c	write(21,*) dsigp00,dsigi00,sigma_eep
350              c	write(21,*) dsigp0,dsigi0,sigma_eep
351              	return
352              
353              	end
354              
355              	subroutine eekeeks(ss,q22,angl,theta,phi,epsi,sigma_eep)
356              
357              	include 'simulate.inc'
358 gaskelld 1.1 
359              	real*8 ss,q22,angl,theta,phi,epsi
360              	real*8 w,skc2,skc,q0,q0c,qr,aflx,aflxl,an,x,sx
361              	real*8 as
362              	double complex z1a,z2a,z3a,z4a,z7a,z8a,ur,ui
363              	real*8 dsigt00,dsigl00,dsigp00,dsigi00
364              	real*8 zf(2,6),sigma_eep
365              	real*8 qv2,qv,qvc
366              c	real*8 ctt
367              	integer pn,pna(3),i
368              
369              c real*4 for compatability with CERNLIB routine fint.
370              	real*4 px(3),pa(50)
371              
372              	real*4 fint
373              
374              	w=sqrt(ss)*1000.
375              	skc2=(w**2-Mk2-targ%Mrec_struck**2)**2-4.*Mk2*targ%Mrec_struck**2
376              	skc2=max(skc2,0.d0)
377              	skc=sqrt(skc2)/2./w
378              	q0=-(-q22-w**2+Mp2)/2./Mp
379 gaskelld 1.1 	q0c=(-q22+q0*Mp)/w
380              	qr=sqrt(q22)/q0c
381              	qv2=q22+q0**2
382              	qv=sqrt(qv2)
383              	qvc=Mp/w*qv
384              C ctt not used, and causes floating exception if skc=0 (i.e. skc2<0)
385              C	ctt=pi/skc/qvc*1.d+06
386              	aflx=skc/2./w/(w**2-Mp2)*(hbarc)**2*10000.
387              	aflxl=aflx*qr**2
388              	an=angl*180./pi
389              	x=cos(angl)
390              	sx=sin(angl)
391              	pn=3
392              	px(1)=ss
393              	px(2)=q22/1.d+06
394              	px(3)=an
395              	pna(1)=20
396              	pna(2)=10
397              	pna(3)=19
398              	as=0.0
399              	pa(1)=2.851
400 gaskelld 1.1 	pa(2)=2.898
401              	pa(3)=2.945
402              	pa(4)=2.991
403              	pa(5)=3.038
404              	pa(6)=3.085
405              	pa(7)=3.132
406              	pa(8)=3.320
407              	pa(9)=3.507
408              	pa(10)=3.695
409              	pa(11)=3.883
410              	pa(12)=4.070
411              	pa(13)=4.258
412              	pa(14)=4.446
413              	pa(15)=4.633
414              	pa(16)=4.821
415              	pa(17)=5.009
416              	pa(18)=5.196
417              	pa(19)=5.384
418              	pa(20)=5.572
419              	pa(21)=0.0
420              	pa(22)=0.250
421 gaskelld 1.1 	pa(23)=0.376
422              	pa(24)=0.520
423              	pa(25)=0.750
424              	pa(26)=1.000
425              	pa(27)=1.250
426              	pa(28)=1.500
427              	pa(29)=1.750
428              	pa(30)=2.000
429              	do i=31,49
430              	  pa(i)=as
431              	  as=as+10.
432              	enddo
433              	zf(1,1)=dble(fint(pn,px,pna,pa,zsrff1))
434              	zf(1,2)=dble(fint(pn,px,pna,pa,zsrff2))
435              	zf(1,3)=dble(fint(pn,px,pna,pa,zsrff3))
436              	zf(1,4)=dble(fint(pn,px,pna,pa,zsrff4))
437              	zf(1,5)=dble(fint(pn,px,pna,pa,zsrff5))
438              	zf(1,6)=dble(fint(pn,px,pna,pa,zsrff6))
439              	zf(2,1)=dble(fint(pn,px,pna,pa,zsiff1))
440              	zf(2,2)=dble(fint(pn,px,pna,pa,zsiff2))
441              	zf(2,3)=dble(fint(pn,px,pna,pa,zsiff3))
442 gaskelld 1.1 	zf(2,4)=dble(fint(pn,px,pna,pa,zsiff4))
443              	zf(2,5)=dble(fint(pn,px,pna,pa,zsiff5))
444              	zf(2,6)=dble(fint(pn,px,pna,pa,zsiff6))
445              
446              	ur=(1.D0,0.D0)
447              	ui=(0.D0,1.D0)  
448              
449              	z1a=zf(1,1)*ur+zf(2,1)*ui
450              	z2a=zf(1,2)*ur+zf(2,2)*ui
451              	z3a=zf(1,3)*ur+zf(2,3)*ui
452              	z4a=zf(1,4)*ur+zf(2,4)*ui
453              	z7a=zf(1,5)*ur+zf(2,5)*ui
454              	z8a=zf(1,6)*ur+zf(2,6)*ui
455              c
456              c	calculate free cross section using Saghai's form.
457              c
458              	dsigt00=aflx*(abs(z1a)**2+abs(z2a)**2
459                   &		+2.*real(conjg(z1a)*z2a)*x+0.5*sx**2*(abs(z3a)**2
460                   &		+abs(z4a)**2+2.*real(conjg(z1a)*z4a-conjg(z2a)*z3a
461                   &		+conjg(z3a)*z4a*x)))
462              	dsigl00=aflxl*epsi*(abs(z7a)**2+abs(z8a)**2
463 gaskelld 1.1      &		+2.*real(conjg(z7a)*z8a)*x)
464              	dsigp00=aflx*epsi*(sin(theta))**2*cos(2.*phi)*
465                   &		(0.5*abs(z3a)**2+0.5*abs(z4a)**2+real(conjg(z1a)*z4a-
466                   &		conjg(z2a)*z3a+conjg(z3a)*z4a*x))
467              	dsigi00=aflx*sqrt(2.*qr**2*epsi*(1.+epsi))*
468                   &		sin(theta)*cos(phi)*real((z7a)*(conjg(z3a)-conjg(z2a)+
469                   &		conjg(z4a)*x)+z8a*(conjg(z1a)+conjg(z3a)*x+conjg(z4a)))
470              
471              	sigma_eep=dsigt00+dsigl00+dsigp00+dsigi00
472              	return
473              
474              	end
475              
476              	subroutine phspwght(delta,theta,phi,weight)
477              
478              	include 'simulate.inc'
479              
480              	real*8 delta,theta,phi
481              	real*8 weight
482              	integer pn,pna(3),i
483              	real*8 ps,ts,ds
484 gaskelld 1.1 
485              c real*4 for compatability with CERNLIB routine fint.
486              	real*4 px(3),pa(78)
487              
488              	real*4 fint
489              
490              	if(electron_arm.eq.1 .and. hadron_arm.eq.2)then	!e- in HMS, K in SOS
491              	  pn=2
492              	  px(1)=theta
493              	  px(2)=phi
494              	  pna(1)=20
495              	  pna(2)=50
496              	  ts=-0.03325
497              	  ps=-0.0735
498              	  do i=1,20
499              	    pa(i)=ts
500              	    ts=ts+0.0035
501              	  enddo
502              	  do i=21,70
503              	    pa(i)=ps
504              	    ps=ps+0.003
505 gaskelld 1.1 	  enddo
506              	  weight=dble(fint(pn,px,pna,pa,weightc))
507              	else if (electron_arm.eq.2 .and. hadron_arm.eq.1) then !e- in SOS,K in HmS
508              	  pn=3
509              	  px(1)=delta
510              	  px(2)=theta
511              	  px(3)=phi
512              	  pna(1)=8
513              	  pna(2)=40
514              	  pna(3)=30
515              	  ts=-0.063375
516              	  ps=-0.0435
517              	  ds=-17.5
518              	  do i=1,8
519              	    pa(i)=ds
520              	    ds=ds+5.
521              	  enddo
522              	  do i=9,48
523              	    pa(i)=ts
524              	    ts=ts+0.00325
525              	  enddo
526 gaskelld 1.1 	  do i=49,78
527              	    pa(i)=ps
528              	    ps=ps+0.003
529              	  enddo
530              	  weight=dble(fint(pn,px,pna,pa,weightd))
531              	else
532              	  write(6,*) 'electron_arm=',electron_arm,' and hadron_arm=',hadron_arm
533              	  write(6,*) 'eekeek.f has a phase space factor that is only defined for'
534              	  write(6,*) 'hms&sos case.  Need to update for other spectrometers.'
535              	  stop
536              	endif
537              	weight=max(weight,0.01D00)
538              	weight=max(100.D00/weight,1.0D00)
539              	return
540              	end