EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
radgen.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file radgen.f
1 
2 C **********************************************************************
3 C
4 C RADATA
5 C
6 C **********************************************************************
7 
8  BLOCK DATA radata
9  implicit real*8(a-h,o-z)
10 
11  include "cmpcom.inc"
12  include "ppicom.inc"
13 
14  data
15  + amm/2.7928456d0/,
16  + amn/-1.913148d0/,
17  + chbar/.197328d0/,
18  + barn/.389379d6/,
19  + aml/.511000d-3/,
20  + aml2/.261112d-6/,
21  + al2/.522240d-6/,
22  + isf20/8/,
23  + pi/3.1415926d0/,
24  + pi2/9.869604d0/,
25  + alfa/.729735d-2/,
26  + amc2/1.151857d0/,
27  + amp/.938272d0/,
28  + amh/.938272d0/,
29  + i2/1,1,1,2,3,3,1,2/,
30  + i1/3,3,4,4,3,3,3,3/
31 
32  end
33 
34 C **********************************************************************
35 C
36 C RADGEN
37 C
38 C **********************************************************************
39 
40  SUBROUTINE radgen(e1,q2l,nul,ys,xs,phil,PhRAD,q2tr,anutr,WEIGHT)
41  IMPLICIT NONE
42 ***********************************************************
43 * INPUT :
44 * e1 : beam energy
45 * q2,nu : uncorrected kinematics of scattered lepton
46 * phil : azimuthal angle of scattered lepton
47 * OUTPUT :
48 * PhRAD(4) : real rad. photon 4 vector
49 * q2tr : true Q2
50 * Utr : true U
51 ***********************************************************
52 
53  include "cmpcom.inc"
54  include "radgen.inc"
55  include "phiout.inc"
56  include "tailcom.inc"
57  include "radgenkeys.inc"
58 
59  REAL phrad(*),q2tr,anutr
60  REAL e1,q2l,nul,phil,xs,ys,weight
61  INTEGER itagen
62 *
63 
64 * reset proton mass (abr)
65 * necessary to get correct kinematics after first elastic event !
66  amp=0.938272d0
67 
68 C... Don't recalculate x and y, but use the ones already calculated
69 C ys=geny
70 C xs=genx
71 C ys=nul/e1
72 C xs=q2l/(2.*amp*nul)
73 
74 C write(*,*)ys, q2l, xs, nul, phil, e1
75 * Do we want to use a lookup table?
76  if(ixytest.ge.0)then
77  if(kill_elas_res.ne.2)then
78  call itafun(ys,xs,plrun,pnrun,ixytest,itagen)
79  else
80  itagen=0
81  sigcor=1.
82  endif
83  else
84  itagen=-1
85  endif
86 
87  if(itagen.ne.0.and.abs(ixytest).ne.2)then
88  iphi=0
89  call mpolrad(e1,ys,xs,1.,plrun,pnrun,itagen)
90  if(ixytest.eq.1)ixytest=-2
91  endif
92  call radgam_pol(e1,ys,xs,phil,ixytest,itagen,q2tr,anutr)
93 
94 * set kinematics of possibly radiated real gamma
95  phrad(1)=sngl(dplabg(1))
96  phrad(2)=sngl(dplabg(2))
97  phrad(3)=sngl(dplabg(3))
98 C to avoid rounding errors !
99 CCC phrad(4)=ys*e1-anutr
100  phrad(4)=sqrt(phrad(1)*phrad(1)+phrad(2)*phrad(2)+
101  + phrad(3)*phrad(3))
102 
103 * set the weight
104  if(kill_elas_res.eq.2)then
105  weight=1.
106  else
107  weight = sigcor
108  endif
109 
110  return
111  end
112 
113 C **********************************************************************
114 C
115 C PHIDST_POL
116 C
117 C **********************************************************************
118 
119  subroutine phidst_pol
120 * -----------------
121 *
122 * calculate phi-distribution of bremsstrahlung photons
123 * needed for correction of hadronic distributions
124 *
125  include "double.inc"
126  include "phiout.inc"
127  include "cmpcom.inc"
128  include "radgen.inc"
129  include "ppicom.inc"
130 C
131  dimension db(3),dx(31),dy(31),dz(31)
132  external dfufi_pol
133  data ndim /31/
134 C
135 C boundaries for 2 integrations
136 C
137  db(1)=0.d0
138  db(2)=pi/10.d0
139  db(3)=pi
140  k=0
141  dsum=0.d0
142  do i=1,2
143  dbnd1=db(i)
144  dbnd2=db(i+1)
145  call radquad_pol(dfufi_pol,dbnd1,dbnd2,dx,dy,dz,dh,ndim)
146  do j=1,ndim
147  if(i.eq.1.or.j.ne.1) k=k+1
148  dphi(k)=dx(j)
149  ddephi(k)=dh
150  dsumph(k)=dsum+dz(j)
151  enddo
152  dsum=dsum+dz(ndim)
153  enddo
154  kmp=k
155  return
156  end
157 
158 C **********************************************************************
159 C
160 C RADGAM_POL
161 C
162 C **********************************************************************
163 
164  subroutine radgam_pol(erad,yrad,xrad,genphi,
165  + ixytest,itagen,q2tr,anutr)
166 
167  include "double.inc"
168  include "phiout.inc"
169  include "radgen.inc"
170  include "cmpcom.inc"
171  include "sxycom.inc"
172  include "tailcom.inc"
173  include "amf2com.inc"
174  include "density.inc"
175 
176  real*8 podinl
177  real temp ! Temp variable for a denom
178 
179 C selection of mass dmj and mass index iadmj
180  r1=rlu(0)
181  if(ixytest.lt.0)then
182  rr1=r1*(tbor+tine+tpro+tnuc)
183 C print *,tine,tpro,tnuc
184  if (rr1.gt.(tine+tpro+tnuc)) then
185  ita=0
186  elseif(rr1.gt.(tpro+tnuc)) then
187  ita=1
188  sicurr=rr1-tpro-tnuc
189  call conk2(dble(erad),dble(xrad),dble(yrad),1)
190  elseif(rr1.gt.tnuc)then
191  ita=3
192  sicurr=rr1-tnuc
193  call conk2(dble(erad),dble(xrad),dble(yrad),1)
194  else
195  ita=2
196  sicurr=rr1
197  call conk2(dble(erad),dble(xrad),dble(yrad),2)
198  endif
199  else
200  xs=dble(xrad)
201  ys=dble(yrad)
202  ita=itagen
203  if(ita.eq.1)sicurr=r1*tine
204  if(ita.eq.2)sicurr=r1*tnuc
205  if(ita.eq.3)sicurr=r1*tpro
206  if(ita.eq.2)then
207  call conk2(dble(erad),dble(xrad),dble(yrad),2)
208  else
209  call conk2(dble(erad),dble(xrad),dble(yrad),1)
210 C print *,an,xs,ys,ita,s,x
211  endif
212  isf1=1
213  isf2=isf20
214  isf3=1
215 
216 C print *,an,xs,ys,erad,xrad,yrad
217 C sicurr=scgen
218  endif
219 
220  if(ita.ne.0)then
221 *
222 *--- selection of theta angle
223 *
224 * maximum value of dsitkm can be smaller than 1 !
225 * ---> use modified random number (as in old generator)
226 * abr (03.01.02)
227 * print *,dsitkm(ndxtkm(ita),ita)
228  r1=r1*dsitkm(ndxtkm(ita),ita)
229  do ktk=1,ndxtkm(ita)
230 C print *,dsitkm(ktk,ita),r1,ndxtkm(ita)
231  if (r1.le.dsitkm(ktk,ita)) then
232  if (ktk.eq.1) then
233  dtk1=r1/dsitkm(ktk,ita)
234  else
235  dtk1=(r1-dsitkm(ktk-1,ita))/
236  + (dsitkm(ktk,ita)-dsitkm(ktk-1,ita))
237  endif
238  dtk=dcvtkm(ktk,ita)+(dtk1-0.5d0)*ddetkm(ktk,ita)
239  goto 30
240  endif
241  enddo
242  dtk=dcvtkm(ndxtkm(ita),ita)
243  30 continue
244 * up to this point identical to old generator
245 * exception: calculation of integrand
246 *
247  taout=(sx-sqly*cos(dtk))/ap2
248 C print *,taout,dtk,ita,ktk,r1
249 C print *,sx,sqly,ap2
250  if(ita.eq.1)then
251  taa=taout
252  iphi=0
253  call tails(taout,tm)
254  rcal=ap*demin
255  rmax=(w2-amc2)/(1.+taout)
256  do krr=1,nrr
257  ddeler(krr,ktk)=(rmax-rcal)/nrr
258  drcurr(krr,ktk)=rcal+ddeler(krr,ktk)*(krr-0.5)
259 C print *,krr,ktk,drcurr(krr,ktk)
260  if(krr.eq.1)then
261  dsigmr(krr,ktk)=podinl(drcurr(krr,ktk))
262  else
263  dsigmr(krr,ktk)=dsigmr(krr-1,ktk)+podinl(drcurr(krr,ktk))
264  endif
265 C print *,dsigmr(krr,ktk),drcurr(krr,ktk),taout
266  enddo
267  r2=rlu(0)
268  dsircu=r2*dsigmr(nrr,ktk)
269  do krr=1,nrr
270  if (dsircu.le.dsigmr(krr,ktk)) then
271  if(krr.eq.1)then
272  drr1=dsircu/dsigmr(krr,ktk)
273  else
274  drr1=(dsircu-dsigmr(krr-1,ktk))
275  + /(dsigmr(krr,ktk)-dsigmr(krr-1,ktk))
276  endif
277  drr=drcurr(krr,ktk)+(drr1-0.5d0)*ddeler(krr,ktk)
278  goto 20
279  endif
280  enddo
281  drr=drcurr(nrr,ktk)
282  20 continue
283  dom=drr/ap
284  else
285  dom=(sx-y)/ap/(1.+taout)
286  endif
287  rrout=ap*dom
288 
289 * the following is again the same as in the old generator (abr 03.01.02)
290 *
291 *---selection of phi angle
292 *
293 *
294 C print *,rrout,taout
295 C pause
296  iphi=1
297  call phidst_pol
298 *
299  r3=rlu(0)
300  dphpoi=r3*dsumph(kmp)
301  do kph=2,kmp
302  if (dphpoi.le.dsumph(kph)) then
303 * avoid division by zero (abr)
304 * dphk1=(dphpoi-dsumph(kph-1))/
305 * + (dsumph(kph)-dsumph(kph-1))
306  temp = dsumph(kph) - dsumph(kph-1)
307  if (temp .eq. 0) temp = 1e-20
308  dphk1 = (dphpoi - dsumph(kph-1)) / temp
309  dphk=dphi(kph)+(dphk1-1.0d0)*ddephi(kph)
310  goto 40
311  endif
312  enddo
313  dphk=dphi(kmp)
314  40 continue
315  r4=rlu(0)
316  if (r4.gt.0.5) dphk=-dphk
317 
318 *
319 * radiative photon
320 *
321  deg=dom
322  dthg=dtk
323  dphig=dphk
324  sigma_total=sngl(tbor+tine+tpro+tnuc)
325 
326  q2tr=y+rrout*taout
327  anutr=sx/ap-dom
328 C write(6,*) 'ita,deg,dthg,dphig,q2tr,anutr,xtr',
329 C & ita,deg,dthg,dphig,q2tr,anutr,
330 C & q2tr/(2.*amp*anutr)
331 
332  dgpz=deg*dcos(dthg)
333  dgpxy=deg*dsin(dthg)
334  dgpx=dgpxy*dcos(dphig)
335  dgpy=dgpxy*dsin(dphig)
336 *
337 *---momentum components in the LAB-system:
338 *---two rotations needed - first within the scattering plane
339 *---around the y-axis and second around the new z-axis (beam
340 *---direction) by Phi of the scattering plane
341 *
342  dgplx=-dgpz*dsts+dgpx*dcts
343  dgply=dgpy
344  dgplz=dgpz*dcts+dgpx*dsts
345  dcphi=dcos(dble(genphi))
346  dsphi=dsin(dble(genphi))
347  dplabg(1)=dgplx*dcphi-dgply*dsphi
348  dplabg(2)=dgplx*dsphi+dgply*dcphi
349  dplabg(3)=dgplz
350 
351  else
352 
353  q2tr=2d0*amh*erad*xrad*yrad
354  anutr=yrad*erad
355 C write(*,'(i5,5f8.3)')ita,xs,ys,q2tr,anutr,sigcor
356 
357  dplabg(1)=0.0d0
358  dplabg(2)=0.0d0
359  dplabg(3)=0.0d0
360 
361  endif
362 
363 * reset kinematics (abr)
364 * call conk2(dble(erad),dble(xrad),dble(yrad),1)
365 * now done at beginning of sr radgen
366 
367  end
368 
369 C **********************************************************************
370 C
371 C DFUFI_POL
372 C
373 C **********************************************************************
374 
375  double precision function dfufi_pol(dx)
376 * -----------------------------------
377 * needed for correction of hadronic distributions
378  implicit real*8(a-h,o-z)
379 
380  include "cmpcom.inc"
381  include "sxycom.inc"
382  include "tailcom.inc"
383  include "amf2com.inc"
384  include "ppicom.inc"
385  include "radgen.inc"
386 
387  phipoi=dx
388  taa=taout
389  call tails(taout,tm)
390  dfufi_pol=an*alfa/pi*podinl(rrout)
391 
392  return
393  end
394 
395 C **********************************************************************
396 C
397 C RADQUAD_POL
398 C
399 C **********************************************************************
400 
401  subroutine radquad_pol(dfunct,dlow,dup,dx,dy,dz,dh,ndim)
402 * -------------------------------------------------
403 *
404  include "double.inc"
405  dimension dx(31),dy(31),dz(31)
406  external dfunct
407 *
408  dsum2=0.d0
409  if (ndim.gt.1) then
410  dh=(dup-dlow)/float(ndim-1)
411  do i=1,ndim
412  daux=dlow+dh*float(i-1)
413  dx(i)=daux
414  dy(i)=dfunct(daux)
415  enddo
416  do i=2,ndim
417  dsum1=dsum2
418  dsum2=dsum2+.5d0*(dx(i)-dx(i-1))*(dy(i)+dy(i-1))
419  dz(i-1)=dsum1
420  enddo
421  dz(ndim)=dsum2
422  elseif (ndim.eq.1) then
423  dz(ndim)=dsum2
424  endif
425 *
426  return
427  end
428 
429 C **********************************************************************
430 C
431 C MPOLRAD
432 C
433 C **********************************************************************
434 
435  subroutine mpolrad(e1curr,yscurr,xscurr
436  + ,uncurr,plcurr,pncurr,itagen)
437  implicit real*8(a-h,o-z)
438 
439  include "cmpcom.inc"
440  include "sxycom.inc"
441  include "ppicom.inc"
442  include "tailcom.inc"
443  include "radgen.inc"
444  include "radgenkeys.inc"
445  include "deltacom.inc"
446  include "density.inc"
447  include "pypars.inc"
448 
449  real e1curr,yscurr,xscurr,uncurr,plcurr,pncurr
450 
451  dimension tai(5),si(2,3),si0(2,3),tls(2,3,4)
452 
453  e1=e1curr
454  xs=xscurr
455  ys=yscurr
456  pl=plcurr
457  pn=-pncurr
458  un=uncurr
459 
460  if (abs(pncurr).eq.2.) then
461 c negative polarization is twice as large as positive one
462  qn=pncurr*1./2.
463  if(qn.lt.0) qn=2.*qn
464  else
465  qn=0.
466  endif
467 
468  call conk2(e1,xs,ys,1)
469 c
470 c
471 c delta is factorizing part of virtual and real leptonic bremsstrahlung
472 c
473  call deltas(delta,delinf,tr)
474 
475  il=1
476  in=1
477 
478  si(il,in)=0d0
479  si0(il,in)=0d0
480  tls(il,in,1)=0d0
481  tls(il,in,2)=0d0
482  tls(il,in,3)=0d0
483  tls(il,in,4)=0d0
484 
485  isf1=1
486  isf2=isf20
487  isf3=1
488 
489  do ii=1,4
490  tai(ii)=0d0
491  enddo
492  sib=0d0
493  sia=0d0
494 
495  if(itagen.eq.-1)then
496  ita1=1
497  ita4=4
498  else
499  ita1=itagen
500  ita4=itagen
501  endif
502 
503  do 30 ita=ita1,ita4
504 c
505 c sib is born cross section with polarized initial
506 c lepton and proton
507 c sia is contribution of anomalous magnetic moment.
508 c
509 * for photoproduction set elastic and quasielastic tails to zero
510 * abr (19.02.04)
511  if (ntx.eq.ntpho .and. (ita.eq.2.or.ita.eq.3)) then
512  tai(2)=0.0d0
513  tai(3)=0.0d0
514  goto 30
515  endif
516 
517 C ECA to use Antjes trick to kill the elastic and quasielastic tails
518 C always if we run pythia with radiative corrections
519 c if (MSTP(199).eq.1) then
520 c tai(2)=0.0d0
521 c tai(3)=0.0d0
522 c goto 30
523 c endif
524 
525  if(ita.eq.2.and.kill_elas_res.eq.1)goto30
526  if(ita.eq.3.and.ire.le.1)then
527  tai(3)=0d0
528 C write(9,'('' tai = 0.0 '')')
529  goto 30
530  end if
531  if(ita.eq.1)then
532  call bornin(sib,sia)
533  endif
534 c
535 c tai(1),tai(2),tai(3) are contributions of radiative tails:
536 c 1 - inelastic
537 c 2 - elastic
538 c 3 - quasielastic
539 c
540  if(ita.eq.2) call conk2(e1,xs,ys,ita)
541 
542  if(kill_elas_res.ne.2) call qqt(sib,tai(ita))
543 
544 * what is this ??? (abr 17.10.01)
545 * shouldnot this if statement be not equal here ???
546 * abr if(ita.eq.2) call conk2(e1,xs,ys,1)
547 * if(ita.ne.2) call conk2(e1,xs,ys,1)
548 * no, purpose is to reset kinematics, I think (abr 03.01.02)
549 * thus the original code is correct
550  if(ita.eq.2) call conk2(e1,xs,ys,1)
551 
552  30 continue
553 C extai1=exp(alfa/pi*delinf)
554 C extai2=((sx-y/rtara)**2/s/(s-y/rtara))**tr
555 C extai3=((sx-y)**2/s/(s-y))**tr
556 
557  extai1=1.
558  extai2=1.
559  extai3=1.
560  delinf=0.
561 
562 C sig=sib*extai1*(1.+alfa/pi*(delta-delinf))+sia
563 C + +tai(4)+tai(5)
564 cilyichev only one-loop without multy-soft photon emission
565 c sig=sib*redfac*(1.+vertex+vac+small)+sia
566  sig=sib*(1.+log(redfac)+vertex+vac+small)+sia
567  + +tai(4)
568  + +tai(1)+(tai(2)*extai2+tai(3)*extai3)/rtara
569 
570  si(il,in)=si(il,in)+sig
571  si0(il,in)=si0(il,in)+sib
572  tls(il,in,1)=tls(il,in,1)+tai(1)
573  tls(il,in,2)=tls(il,in,2)+tai(2)*extai2/rtara
574  tls(il,in,3)=tls(il,in,3)+tai(3)*extai3/rtara
575 C tls(il,in,4)=sib*extai1*(1.+alfa/pi*(delta-delinf))+sia
576 C + +tai(4)+tai(5)
577  tls(il,in,4)=sib*redfac*(1.+vertex+vac+small)+sia
578  + +tai(4)
579 
580 C write(*,'(1x,8g11.4)')xs,ys,si0(il,in),si(il,in)
581 C + ,tls(il,in,1),tls(il,in,2),tls(il,in,3)
582 C + ,tls(il,in,4)
583  if(itagen.eq.-1.or.itagen.eq.1)then
584 C if(abs(tine-tls(il,in,1))/tine.gt.5d-2)
585 C + write(*,*)' tine',tine,tls(il,in,1),itagen
586  tine=tls(il,in,1)
587  endif
588  if(itagen.eq.-1.or.itagen.eq.2)then
589 C if(abs(tnuc-tls(il,in,2))/tnuc.gt.5d-2)
590 C + write(*,*)' tnuc',tnuc,tls(il,in,2),itagen
591  tnuc=tls(il,in,2)
592  endif
593  if(itagen.eq.-1.or.itagen.eq.3)then
594 C if(abs(tpro-tls(il,in,3))/tpro.gt.5d-2)
595 C + write(*,*)' tpro',tnuc,tls(il,in,3),itagen
596  tpro=tls(il,in,3)
597  endif
598  if(itagen.eq.-1)then
599  tbor=tls(il,in,4)
600  sigrad=si(il,in)
601  sig1g=si0(il,in)
602  endif
603 C #ifdef TEST
604 C write(*,'(1x,6g11.4)')sig1g,sigrad,tine,tnuc,tpro,tbor
605 C #endif
606  sigcor=sigrad/sig1g
607  if(kill_elas_res.eq.2)sigcor=1.
608 
609 C #ifdef TEST
610 C print *,' kill_elas_res=',kill_elas_res,amc2
611 C print *,' sig1g=',sig1g,sib
612 C print *,' sigrad=',sigrad,sig
613 C print *,' tine =',tine
614 C print *,' tai(4) =',tai(4)
615 C print *,' tnuc =',tnuc
616 C print *,' tpro =',tpro
617 C print *,' tbor =',tbor
618 C print *,' sig1g=',sig1g
619 C print *,' sigcor=',sigcor
620 C print *,' delta=',alfa/pi*delta
621 C print *,' vac =',vac
622 C print *,' vertex=',vertex
623 C print *,' small=',small
624 C print *,' tai(4)= ',sib*log(redfac)
625 C sigt1=sib*alfa/pi*delta+tai(5)
626 C sigt2=sib*(log(redfac)+vertex+vac+small)
627 C sigt3=sib*alfa/pi*(delta+delta5)
628 C print *,' sigt=',sigt1,sigt2,sigt3
629 C #endif
630 * print *,' kill_elas_res=',kill_elas_res,amc2
631 * print *,' sig1g=',sig1g,sib
632 * print *,' sigrad=',sigrad,sig
633 * print *,' tine =',tine
634 * print *,' tai(4) =',tai(4)
635 * print *,' tnuc =',tnuc
636 * print *,' tpro =',tpro
637 * print *,' tbor =',tbor
638 * print *,' sig1g=',sig1g
639 * print *,' sigcor=',sigcor
640 * print *,' delta=',alfa/pi*delta
641 * print *,' vac =',vac
642 * print *,' vertex=',vertex
643 * print *,' small=',small
644 * print *,' redfac= ',redfac
645 
646  end
647 
648 C **********************************************************************
649 C
650 C CONK2
651 C
652 C **********************************************************************
653 
654  subroutine conk2(e1,xs,ys,iittaa)
655  implicit real*8(a-h,o-z)
656 
657  include "cmpcom.inc"
658 
659  if(iittaa.eq.2)then
660  amp=amt
661  else
662  amp=amh
663  endif
664 C print *,amp,amh,iittaa
665  call conkin(e1,xs*amh/amp,ys)
666 
667  return
668  end
669 
670 C **********************************************************************
671 C
672 C CONKIN
673 C
674 C **********************************************************************
675 
676  subroutine conkin(e1,xss,yss)
677 
678 c set of kinematical constants
679  implicit real*8(a-h,o-z)
680 
681  include "cmpcom.inc"
682  include "polcom.inc"
683  include "sxycom.inc"
684  include "ppicom.inc"
685  include "radgen.inc"
686 cilyichev additional include
687  include "mc_set.inc"
688 
689  ap=2.*amp
690  amp2=amp**2
691  ap2=2.*amp**2
692  s=ap*e1
693  x=s*(1-yss)
694  sx=s-x
695  sxp=s+x
696  y=s*xss*yss
697  ym=y+al2
698  tpl=s**2+x**2
699  tmi=s**2-x**2
700  w2=amp2+s-y-x
701  als=s*s-al2*ap2
702  alx=x*x-al2*ap2
703  alm=y*y+4.*aml2*y
704  aly=sx**2+4.*amp2*y
705  sqls=dsqrt(als)
706  sqlx=dsqrt(alx)
707  sqly=dsqrt(aly)
708  sqlm=dsqrt(alm)
709  allm=dlog((sqlm+y)/(sqlm-y))/sqlm
710 
711  coe=xs/e1/1d3
712 
713  axy=pi*(s-x) * coe
714  an=2.*alfa**2/sqls*axy*barn*amh/amp
715  tamin=(sx-sqly)/ap2
716  tamax=(sx+sqly)/ap2
717 
718  dcts=(s*sx + ap2*y)/sqly/sqls
719 * avoid values above 1.0 (for x values below 1e-07, abr 18.2.04)
720 * dsts=sin( acos(dcts) )
721  dsts=sin( acos(min(max(dcts,-1.0),+1.0)) )
722 
723  as=s/2./aml/sqls
724  bs=0.
725  cs=-aml/sqls
726 cilyichev instead of
727 c ae=amp/sqls
728 c be=0.
729 c ce=-s/ap/sqls
730 c we need
731  if ((mcset_ptarget(1:1).eq.'L').or.
732  + (mcset_ptarget(1:2).eq.'DT')) then
733  ae=amp/sqls
734  be=0.
735  ce=-s/ap/sqls
736  elseif(mcset_ptarget(1:1).eq.'T')then
737  sqn=dsqrt(s*x*y-aly*aml2-amp2*y*y)
738  ae=(-s*x+ap2*ym)/sqls/sqn/2.
739  be=sqls/sqn/2.
740  ce=-(s*y+al2*sx)/sqls/sqn/2.
741  endif
742 
743  apq=-y*(ae-be)+ce*sx
744  apn=(y+4.*aml2)*(ae+be)+ce*sxp
745  dk2ks=as*ym+al2*bs+cs*x
746  dksp1=as*s+bs*x+cs*ap2
747  dapks=2.*(al2*(as*ae+bs*be)+ap2*cs*ce+ym*(as*be+bs*ae)
748  + +s*(as*ce+cs*ae)+x*(bs*ce+cs*be))
749 
750  return
751  end
752 
753 C **********************************************************************
754 C
755 C BORNIN
756 C
757 C **********************************************************************
758 
759  subroutine bornin(sibor,siamm)
760 c
761 c sibor is born cross section with polarized initial
762 c lepton and polarized target
763 c siamm is contribution of anomalous magnetic moment.
764 C the cross section is calculated as dsigma/dlogQ2dnu to fit
765 C to the HERMES-disng generation
766 
767 c
768  implicit real*8(a-h,o-z)
769 
770  include "cmpcom.inc"
771  include "polcom.inc"
772  include "sxycom.inc"
773  include "ppicom.inc"
774  include "tailcom.inc"
775 
776  dimension sfm0(8),tm(8)
777 
778  call strf(0d0,0d0,sfm0)
779  tm(1)=-(2.*aml2-y)
780  tm(2)=(-(amp2*y-s*x))/(2.*amp2)
781  tm(3)=(2.*(apq*dk2ks-dapks*y)*aml)/amp
782  tm(4)=apq/amp*(-(dk2ks*sx-2.*dksp1*y)*aml)/amp2
783  tm(7)=(-(4.*aml2+3.*apn**2-3.*apq**2+y))/2.
784  tm(8)=apq/amp*(-3.*(apn*sxp-apq*sx))/(2.*ap)
785  ek=(3.*apq**2-y)/amp2
786  tm(5)=-ek*tm(1)
787  tm(6)=-ek*tm(2)
788  ssum=0.
789  do 1 isf=isf1,isf2,isf3
790  ppol=1.
791  if(isf.eq.3.or.isf.eq.4)ppol=pl*pn
792 C if(isf.eq.3.or.isf.eq.4)ppol=-pn
793  if(isf.ge.5)ppol=qn/6
794  ssum=ssum+tm(isf)*sfm0(isf)*ppol
795  1 continue
796  sibor=ssum*an/y/y*2.
797 c
798 c formula (4) of kukhto and shumeiko paper
799 c
800 cc res1=amp*ww1*(y+4.*aml2)-ww2*(s+x)**2/4./amp
801 cc siamm=alfa/pi*al2*allm*(sibor+an*res1/y**2)
802  siamm=0.
803 
804  return
805  end
806 
807 C **********************************************************************
808 C
809 C DELTAS
810 C
811 C **********************************************************************
812 
813  subroutine deltas(delta,delinf,tr)
814 c
815 c delta is factorizing part of virtual and real leptonic bremsstrahlung
816 c
817  implicit real*8(a-h,o-z)
818 
819  include "cmpcom.inc"
820  include "sxycom.inc"
821  include "ppicom.inc"
822  include "radgen.inc"
823  include "deltacom.inc"
824 
825 c
826 c am2 : squared masses of charge leptons
827 c
828  dimension am2(3)
829  data am2/.26110d-6,.111637d-1,3.18301d0/
830 
831  del1=-ym*(alm*allm**2/2.+2.*fspen(2d0*sqlm/(y+sqlm))-pi2/2.)/sqlm
832  del2=(3.*y/2.+4.*aml2)*allm-2.
833 
834  suml=0.
835  do 10 i=1,3
836  a2=2.*am2(i)
837  sqlmi=dsqrt(y*y+2.*a2*y)
838  allmi=dlog((sqlmi+y)/(sqlmi-y))/sqlmi
839  suml=suml+2.*(y+a2)*allmi/3.-10./9.+4.*a2*(1.-a2*allmi)/3./y
840  10 continue
841  if(y.lt.1.d0)then
842  aaa = -1.345d-9
843  bbb = -2.302d-3
844  ccc = 4.091
845  elseif(y.lt.64d0)then
846  aaa = -1.512d-3
847  bbb = -2.822d-3
848  ccc = 1.218
849  else
850  aaa = -1.1344d-3
851  bbb = -3.0680d-3
852  ccc = 9.9992d-1
853  endif
854  sumh = -(aaa+bbb*log(1.+ccc*y)) *2*pi/alfa
855  sum=suml+sumh
856 C print *,' vacl_my_old=',alfa/pi*suml
857 C print *,' vach_my_old=',alfa/pi*sumh
858 C print *,' vac_my_old=',alfa/pi*sum
859  sum=vacpol(y)
860 C print *,' vac_my_new=',alfa/pi*sum
861 
862  aj0=2.*(ym*allm-1.)
863 c...elke okay the dlog((w2-amc2) can become negative if w2 is smaller than
864 C amc2/1.151857d0/
865 c write(*,*)w2,aml,amc2,aj0
866  deltai=aj0*dlog((w2-amc2)/aml/dsqrt(w2))
867 c write(*,*)"deltai=",deltai
868 
869  ss=x+y
870  xx=s-y
871  alss=ss**2-2.*w2*al2
872  alxx=xx**2-2.*w2*al2
873  sqlss=dsqrt(alss)
874  sqlxx=dsqrt(alxx)
875  allss=dlog((sqlss+ss)/(-sqlss+ss))/sqlss
876  allxx=dlog((sqlxx+xx)/(-sqlxx+xx))/sqlxx
877  dlm=dlog(y/aml2)
878  sfpr=dlm**2/2.-dlm*dlog(ss*xx/aml2/w2)
879  + -(dlog(ss/xx))**2/2.+fspen((s*x-y*amp2)/ss/xx)-pi2/3.
880  delta0=(ss*allss+xx*allxx)/2.+sfpr
881  delta=deltai+delta0+del1+del2+sum
882  delinf=(dlm-1.)*dlog((w2-amc2)**2/ss/xx)
883  tr=alfa/pi*(dlm-1.)
884 
885  vac=alfa/pi*sum
886  vertex=alfa/pi*del2
887  small_old=alfa/pi*(pi2/6.-fspen(1.-amp2*y/s/x)
888  + + fspen(1.-s/x) + fspen(1.-x/s))
889  small=alfa/pi*(-pi2/6.-log(ss*xx/s/x)*log(sx/y)
890  + +log(ss/s)*log(xx/x)
891  + -1.5*log(s/x)**2+2.*log(xx/x)*log(s/y)+2.*log(ss/s)*log(x/y)
892  + -2.*fspen(-y/x)-2.*fspen(y/s)+2.*fspen(y/sx)+fspen(s*x/ss/xx)
893  + -fspen(s*y/ss/sx)-fspen(x*y/sx/xx))
894 C print *,' small_old=',small_old
895 C print *,' small_new=',small
896  redfac=exp(-alfa/pi*(dlm-1.)*log(s*x/(4.*amp2*demin**2) ))
897  delta5=-(dlm-1.)*log((w2-amc2)**2*s*x/(4.*amp2*demin**2*ss*xx))
898  + -log(xx/x)*log(amp2*y/s**2)-log(ss/s)*log(amp2*y/x**2)
899  + -2.*fspen(-y/x)-2.*fspen(y/s)+2.*fspen(-tamin)
900  + +2.*fspen(-tamax)
901  + -fspen((-y-s*tamin)/xx)
902  + -fspen((-y-s*tamax)/xx)
903  + -fspen(( y-x*tamin)/ss)
904  + -fspen(( y-x*tamax)/ss)
905 
906  return
907  end
908 
909 C **********************************************************************
910 C
911 C VACPOL
912 C
913 C **********************************************************************
914 
915  double precision function vacpol(y)
916 
917  implicit real*8(a-h,o-z)
918 
919  include "ppicom.inc"
920 
921 c
922 c am2 : squared masses of charge leptons
923 c
924  dimension am2(3),a(5),b(5),c(5)
925  data am2/.26110d-6,.111637d-1,3.18301d0/
926  data a/0d0,0d0,0d0,1.2227d-3,1.64178d-3/
927  data b/2.2877d-3,2.51507d-3,2.79328d-3,2.96694d-3,2.92051d-3/
928  data c/4.08041425d0,3.09624477d0,2.07463133d0,1d0,1d0/
929 
930  suml=0.
931  do 10 i=1,3
932  a2=2.*am2(i)
933  sqlmi=dsqrt(y*y+2.*a2*y)
934  allmi=dlog((sqlmi+y)/(sqlmi-y))/sqlmi
935  suml=suml+2.*(y+a2)*allmi/3.-10./9.+4.*a2*(1.-a2*allmi)/3./y
936  10 continue
937 
938  if(y .lt. 4d0)then
939  k=1
940  elseif(y .lt. 16d0)then
941  k=2
942  elseif(y .lt. 100d0)then
943  k=3
944  elseif(y .lt. 8317.44d0)then
945  k=4
946  elseif(y .ge. 8317.44d0)then
947  k=5
948  else
949  stop ' Y<0 in VACPOL'
950  endif
951 
952  sumh = (a(k)+b(k)*log(1.+c(k)*y)) *2.*pi/alfa
953  vacpol=suml+sumh
954 
955  end
956 
957 C **********************************************************************
958 C
959 C FSPENS
960 C
961 C **********************************************************************
962 
963  double precision function fspens(x)
964 c
965 c spence function
966 c
967  implicit real*8(a-h,o-z)
968 
969  f=0.d0
970  a=1.d0
971  an=0.d0
972  tch=1.d-16
973  1 an=an+1.d0
974  a=a*x
975  b=a/an**2
976  f=f+b
977  if(b-tch)2,2,1
978  2 fspens=f
979 
980  return
981  end
982 
983 C **********************************************************************
984 C
985 C FSPEN
986 C
987 C **********************************************************************
988 
989  double precision function fspen(x)
990 
991  implicit real*8(a-h,o-z)
992 
993  data f1/1.644934d0/
994 
995  if(x)8,1,1
996  1 if(x-.5d0)2,2,3
997  2 fspen=fspens(x)
998  return
999  3 if(x-1d0)4,4,5
1000  4 fspen=f1-dlog(x)*dlog(1d0-x+1d-10)-fspens(1d0-x)
1001  return
1002  5 if(x-2d0)6,6,7
1003  6 fspen=f1-.5*dlog(x)*dlog((x-1d0)**2/x)+fspens(1d0-1d0/x)
1004  return
1005  7 fspen=2d0*f1-.5d0*dlog(x)**2-fspens(1d0/x)
1006  return
1007  8 if(x+1d0)10,9,9
1008  9 fspen=-.5d0*dlog(1d0-x)**2-fspens(x/(x-1d0))
1009  return
1010  10 fspen=-.5*dlog(1.-x)*dlog(x**2/(1d0-x))-f1+fspens(1d0/(1d0-x))
1011  return
1012 
1013  end
1014 
1015 C **********************************************************************
1016 C
1017 C QQT
1018 C
1019 C **********************************************************************
1020 
1021  subroutine qqt(bo,tai)
1022 
1023  implicit real*8(a-h,o-z)
1024 
1025  include "cmpcom.inc"
1026  include "sxycom.inc"
1027  include "ppicom.inc"
1028  include "tailcom.inc"
1029  include "radgen.inc"
1030 
1031  external rv2
1032  dimension dbtk(8)
1033  data ep/1d-8/
1034 
1035  dsumtk=0.d0
1036  derrtk=0.d0
1037  isumtk=0
1038  if(ita.eq.1 .or. ita.eq.5)then
1039  tade=(w2-amc2)/(ap*demin) -1.d0
1040  costkm=(sx-ap2*tade)/sqly
1041  if(costkm.ge.1d0)then
1042  tai=0.
1043  return
1044  endif
1045  dtkmax=acos( max(-1d0,costkm ))
1046  else
1047  dtkmax=pi
1048  endif
1049  call intbtk2(dbtk,nbtk,dtkmax)
1050 c integrate each bin by ntk subbins
1051 c ntk=number of bins within the big bin dbtk(i)...dbtk(i+1)
1052  do 10 itk=1,nbtk
1053  call inttk2(isumtk,dbtk(itk),dbtk(itk+1),dsumtk,derrtk)
1054 C write(*,*)itk,isumtk,dsumtk
1055 
1056  10 continue
1057  if(ita.le.3)ndxtkm(ita)=isumtk
1058  tai=dsumtk
1059 
1060  end
1061 
1062 C **********************************************************************
1063 C
1064 C TAILS
1065 C
1066 C **********************************************************************
1067 
1068  subroutine tails(ta,tm)
1069  implicit real*8(a-h,o-z)
1070 
1071  include "cmpcom.inc"
1072  include "polcom.inc"
1073  include "sxycom.inc"
1074  include "ppicom.inc"
1075  include "radgen.inc"
1076  include "bseocom.inc"
1077 
1078  dimension tm(8,6),ajm2(2),ajm3(3),ii(8)
1079  data ii/1,2,3,4,1,2,5,6/
1080 
1081  if(iphi.eq.0)then
1082  b2=(-aly*ta+sxp*sx*ta+2.*sxp*y)/2.
1083  b1=(-aly*ta-sxp*sx*ta-2.*sxp*y)/2.
1084  c1=-(4.*(amp2*ta**2-sx*ta-y)*aml2-(s*ta+y)**2)
1085  c2=-(4.*(amp2*ta**2-sx*ta-y)*aml2-(ta*x-y)**2)
1086  bb=1./sqly
1087  sc1=dsqrt(c1)
1088  sc2=dsqrt(c2)
1089  bi12=(sxp*(sx*ta+2.*y))/(sc1*sc2*(sc1+sc2))
1090  bi1pi2=1./sc2+1./sc1
1091  bis=-b1/sc1/c1+b2/sc2/c2
1092  bir=b2/sc2/c2+b1/sc1/c1
1093  b1i=-b1/aly/sqly
1094  b11i=(3.*b1**2-aly*c1)/2./aly**2/sqly
1095  else
1096 C write(*,*) ta,tamin,tamax,s,x,y,aml2,aly
1097 C...elke: to avoid problems with the sqrt from rounding at the very low Q2
1098 C if using pythia6 ----> y becomes very small
1099  sqrtmb=((ta-tamin)*(tamax-ta)*(s*x*y-(y**2*amp2)-(aml2*aly)))
1100  if (sqrtmb.ge.0) then
1101  sqrtmb=sqrt(sqrtmb)
1102  else
1103  write(*,*)'radgen: ta,tamin,tamax,s,x,y,aml2,aly',
1104  & ta,tamin,tamax,s,x,y,aml2,aly
1105  write(*,*) "radgen: sqrtmb=",sqrtmb
1106  sqrtmb=0
1107  endif
1108 
1109  z1=(y*sxp+ta*(s*sx+ap2*y)-ap*cos(phipoi)*sqrtmb)/aly
1110  z2=(y*sxp+ta*(x*sx-ap2*y)-ap*cos(phipoi)*sqrtmb)/aly
1111  bb=1./sqly/pi
1112  bi12=bb/(z1*z2)
1113  bi1pi2=bb/z2+bb/z1
1114  bis=bb/z2**2+bb/z1**2
1115  bir=bb/z2**2-bb/z1**2
1116  b1i=bb*z1
1117  b11i=bb*z1**2
1118  endif
1119  sps=as+bs
1120  spe=ae+be
1121  ccpe=(ae-be)*ta+2.*ce
1122  ccps=(as-bs)*ta+2.*cs
1123  sis=(2.*bi1pi2*sps+bir*sps*ta+bis*ccps)/2.
1124  sir=( (2.*bi12*sps*ta+bir*ccps+bis*sps*ta))/2.
1125  si12=(bi12*ccps+bi1pi2*sps)/2.
1126  eis=(2.*bi1pi2*spe+bir*spe*ta+bis*ccpe)/2.
1127  eir=( (2.*bi12*spe*ta+bir*ccpe+bis*spe*ta))/2.
1128  ei12=(bi12*ccpe+bi1pi2*spe)/2.
1129  ois=((2.*bi1pi2+bir*ta)*(ccpe*sps+ccps*spe)+(ccpe*ccps+
1130  + spe*sps*ta**2)*bis+8.*bb*spe*sps+4.*bi12*spe*sps*ta**2)/4.
1131  oir=( ((2.*bi12+bis)*(ccpe*sps+ccps*spe)*ta+(ccpe*ccps+
1132  + spe*sps*ta**2)*bir+4.*bi1pi2*spe*sps*ta))/4.
1133  oi12=((ccpe*ccps+spe*sps*ta**2)*bi12+(ccpe*sps+ccps*spe)*
1134  + bi1pi2+4.*bb*spe*sps)/4.
1135  eeis=((ccpe**2+spe**2*ta**2)*bis+8.*bb*spe**2+4.*bi12*spe
1136  + **2*ta**2+4.*bi1pi2*ccpe*spe+2.*bir*ccpe*spe*ta)/4.
1137  eeir=( ((ccpe**2+spe**2*ta**2)*bir+4.*bi12*ccpe*spe*ta+4.
1138  + *bi1pi2*spe**2*ta+2.*bis*ccpe*spe*ta))/4.
1139  eei12=((ccpe**2+spe**2*ta**2)*bi12+4.*bb*spe**2+2.*bi1pi2
1140  + *ccpe*spe)/4.
1141  ei1pi2=(4.*bb*spe+bi12*spe*ta**2+bi1pi2*ccpe)/2.
1142  eei1i2=((ccpe**2+spe**2*ta**2)*bi1pi2+4.*(2.*ccpe-spe*ta)
1143  + *bb*spe+8.*b1i*spe**2+2.*bi12*ccpe*spe*ta**2)/4.
1144  eb=((ccpe-spe*ta)*bb+2.*b1i*spe)/2.
1145  eeb=((ccpe-spe*ta)**2*bb+4.*(ccpe-spe*ta)*b1i*spe+4.*b11i
1146  + *spe**2)/4.
1147  call ffu(1,bb,bis,bir,bi12,bi1pi2,sir,sis,si12
1148  + ,eis,eir,ei12,ei1pi2,ta)
1149  call ffu(2,eb,eis,eir,ei12,ei1pi2,oir,ois,oi12
1150  + ,eeis,eeir,eei12,eei1i2,ta)
1151  call ffu(3,eeb,eeis,eeir,eei12,eei1i2,0d0,0d0,0d0
1152  + ,0d0,0d0,0d0,0d0,ta)
1153  ajm2(1)=apq/amp
1154  ajm2(2)=-1./amp
1155  ajm3(1)=(y-3.*apq**2)/amp2
1156  ajm3(2)=6.*apq/amp2
1157  ajm3(3)=-3./amp2
1158  do 15 i=1,8
1159  do 13 l=1,6
1160  tm(i,l)=0
1161  13 enddo
1162  do 10 k=1,i2(i)
1163  ajk=1.
1164  if(i.eq.4.or.i.eq.8)ajk=ajm2(k)
1165  if(i.eq.5.or.i.eq.6)ajk=ajm3(k)
1166  do 11 j=k,i1(i)+k-1
1167  tm(i,j)=tm(i,j)+tm3(ii(i),j-k+1,k)*ajk
1168  if((i.eq.5.or.i.eq.6).and.k.eq.2)
1169  + tm(i,j)=tm(i,j)+tm3(ii(i),j-k+1,1)*ta/amp2
1170  11 continue
1171  10 continue
1172  15 continue
1173 
1174  return
1175  end
1176 
1177 C **********************************************************************
1178 C
1179 C FFU
1180 C
1181 C **********************************************************************
1182 
1183  subroutine ffu(n,bb,bis,bir,bi12,bi1pi2,sir,sis,si12
1184  + ,eis,eir,ei12,ei1pi2,ta)
1185  implicit real*8(a-h,o-z)
1186 
1187  include "cmpcom.inc"
1188  include "polcom.inc"
1189  include "sxycom.inc"
1190  include "ppicom.inc"
1191  include "bseocom.inc"
1192 
1193  hi2=aml2*bis-ym*bi12
1194  shi2=aml2*sis-ym*si12
1195  ehi2=aml2*eis-ym*ei12
1196  ohi2=aml2*ois-ym*oi12
1197  goto(10,20,30)n
1198  10 continue
1199  tm3(3,1,n)=(8.*(apq*dk2ks-dapks*y)*aml*hi2)/amp
1200  tm3(3,2,n)=(-2.*((2.*(bi12*dk2ks*ta-2.*shi2)*apq+(2.*shi2-
1201  + sir*y+sis*ym)*apn+4.*dapks*hi2*ta)-4.*((2.*ei12-eis)*
1202  + dk2ks-(si12-sis)*apn)*aml2)*aml)/amp
1203  tm3(3,3,n)=(2.*(((2.*si12+sir-sis)*apn*ta-2.*dk2ks*ei12*ta
1204  + -6.*ohi2-oir*y+ois*ym)-4.*aml2*oi12)*aml)/amp
1205  tm3(3,4,n)=(2.*(2.*oi12-oir+ois)*aml*ta)/amp
1206  tm3(5,1,n)=-2.*(4.*aml2+3.*apn**2-3.*apq**2+y)*hi2
1207  tm3(5,2,n)=-2.*(6.*aml2*apn*eir-3.*apn**2*bi12*ta+3.*apn*
1208  + apq*bi1pi2+6.*apq*ehi2+hi2*ta)
1209  tm3(5,3,n)=-(24.*aml2*eei12-6.*apn*ei1pi2-6.*apq*ei12*ta-
1210  + 2.*bb-bi12*ta**2)
1211  20 continue
1212  tm3(4,1,n)=(-4.*(dk2ks*sx-2.*dksp1*y)*aml*hi2)/amp2
1213  tm3(4,2,n)=(((2.*(sxp-2.*sx)*shi2+2.*bi12*dk2ks*sx*ta+8.*
1214  + dksp1*hi2*ta-sir*sxp*y+sis*sxp*ym)-4.*(2.*bi12*dk2ks-bis*
1215  + dk2ks-si12*sxp+sis*sxp)*aml2)*aml)/amp2
1216  tm3(4,3,n)=((((sxp*ta-ym)*sis-(sxp*ta-y)*sir+2.*bi12*dk2ks
1217  + *ta+6.*shi2-2.*si12*sxp*ta)+4.*aml2*si12)*aml)/amp2
1218  tm3(4,4,n)=(-(2.*si12-sir+sis)*aml*ta)/amp2
1219  tm3(6,1,n)=(-3.*(apn*sxp-apq*sx)*hi2)/amp
1220  tm3(6,2,n)=(-3.*(2.*(apn*bir+eir*sxp)*aml2-(2.*bi12*sxp*ta
1221  + -bi1pi2*sx)*apn+(bi1pi2*sxp+2.*hi2)*apq+2.*ehi2*sx))/(2.*
1222  + amp)
1223  tm3(6,3,n)=(-3.*(8.*aml2*ei12-apn*bi1pi2-apq*bi12*ta-ei12*
1224  + sx*ta-ei1pi2*sxp))/(2.*amp)
1225  30 continue
1226  tm3(1,1,n)=-4.*(2.*aml2-y)*hi2
1227  tm3(1,2,n)=4.*hi2*ta
1228  tm3(1,3,n)=-2.*(2.*bb+bi12*ta**2)
1229  tm3(2,1,n)=(((sxp**2-sx**2)-4.*amp2*y)*hi2)/(2.*amp2)
1230  tm3(2,2,n)=(2.*aml2*bir*sxp-4.*amp2*hi2*ta-bi12*sxp**2*ta+
1231  + bi1pi2*sxp*sx+2.*hi2*sx)/(2.*amp2)
1232  tm3(2,3,n)=(2.*(2.*bb+bi12*ta**2)*amp2+4.*aml2*bi12-bi12*
1233  + sx*ta-bi1pi2*sxp)/(2.*amp2)
1234 
1235  return
1236  end
1237 
1238 C **********************************************************************
1239 C
1240 C PODINL
1241 C
1242 C **********************************************************************
1243 
1244  double precision function podinl(r)
1245 c
1246 c integrand (over r )
1247 c
1248  implicit real*8(a-h,o-z)
1249 
1250  include "cmpcom.inc"
1251  include "sxycom.inc"
1252  include "tailcom.inc"
1253  include "ppicom.inc"
1254  include "amf2com.inc"
1255 
1256  dimension sfm(8),iel(8)
1257  data iel/0,2,1,2,2,4,2,3/
1258 
1259  if(ita.ne.5) call strf(taa ,r,sfm)
1260  podinl=0.
1261  do 11 isf=isf1,isf2,isf3
1262  ppol=1.
1263  if(isf.eq.3.or.isf.eq.4)ppol=pl*pn
1264 C if(isf.eq.3.or.isf.eq.4)ppol=-pn
1265  if(isf.ge.5)ppol=qn/6
1266  if(ita.eq.2)ppol=ppol*(amt/amp)**iel(isf)
1267  irrend=i1(isf)+i2(isf)-1
1268  if(ita.eq.5)irrend=1
1269  do 1 irr=1,irrend
1270  if(ita.eq.5)then
1271  pres=-sfm0(isf)/2./y**2/r
1272  else
1273  pp=sfm(isf)
1274  if(irr.eq.1.and.ita.eq.4)pp=pp-sfm0(isf)*(1.+r*taa/y)**2
1275  pres=pp*r**(irr-2)/(y+r*taa)**2/2.
1276  endif
1277  podinl=podinl-tm(isf,irr)*pres*ppol
1278 C print *,' isf=',isf,' irr=',irr,podinl
1279 C write(*,*)' isf=',isf,' irr=',irr,podinl,sfm(isf)
1280  1 continue
1281  11 continue
1282  podinl=podinl*(1.+alfa/pi*vacpol(y+r*taa))
1283 
1284  return
1285  end
1286 
1287 C **********************************************************************
1288 C
1289 C RV2
1290 C
1291 C **********************************************************************
1292 
1293  double precision function rv2(ta)
1294 c
1295 c integrand (over ta )
1296 c
1297  implicit real*8(a-h,o-z)
1298  external podinl
1299 
1300  include "cmpcom.inc"
1301  include "sxycom.inc"
1302  include "tailcom.inc"
1303  include "amf2com.inc"
1304  include "ppicom.inc"
1305  include "radgen.inc"
1306 
1307  taa=ta
1308  call strf(0d0,0d0,sfm0)
1309  call tails(ta,tm)
1310  rmin=1d-8
1311  rcal=ap*demin
1312  rmax=(w2-amc2)/(1.+ta)
1313  if(ita.eq.1)then
1314 C call dqn32(rmin,rmax,podinl,res)
1315  dsumtk=0.d0
1316  derrtk=0.d0
1317 
1318  dbmj2=rmax
1319  dbmj1=rcal
1320  nmj=nrr
1321  dr=(rmax-rcal)/nrr
1322  rcur=rcal-dr*.5d0
1323  dsum1=0.d0
1324 c loop over all bins
1325  do irr=1,nrr
1326  rcur=rcur+dr
1327  drcurr(irr,itkcur) = rcur
1328  ddeler(irr,itkcur) = dr
1329  dsum1=dsum1+podinl(rcur)
1330  dsigmr(irr,itkcur) = dr*dsum1
1331  enddo
1332 
1333  res=dr*dsum1
1334 
1335  elseif(ita.eq.2 .or. ita.eq.3)then
1336  aa=amt/amp
1337  if(ita.eq.3)aa=1.
1338  res=podinl((sx-y/aa)/(1d0+ta/aa))/(1.+ta/aa) /aa**2
1339  elseif(ita.eq.4)then
1340  rend=min(rcal,rmax)
1341  call dqn32(rmin,rend,podinl,res)
1342 C call qunc8(podinl,rmin,rend,1d-4,1d-9,res,er,nn,fl,3500)
1343  elseif(ita.eq.5)then
1344  res=podinl(1.d0/log(rmax/rcal))
1345  endif
1346  rv2=res
1347 C write(9,*)res,dr,rcur
1348 
1349  return
1350  end
1351 
1352 C **********************************************************************
1353 C
1354 C STRF
1355 C
1356 C **********************************************************************
1357 
1358  subroutine strf(ta,rr,sfm)
1359 c
1360 c the programm calculates deep inelastic (ita=1),
1361 c elastic (ita=2), quasielastic (ita=3) structure functions
1362 c in kinematical point (ta,rr).
1363 c rr=sx-tt,
1364 c ta=(t-y)/rr,
1365 c where tt=t+amf2-amp2, amf2 is invarint mass of final hadrons
1366 c
1367  implicit real*8(a-h,o-z)
1368 
1369  include "cmpcom.inc"
1370  include "sxycom.inc"
1371  include "tailcom.inc"
1372  include "radgenkeys.inc"
1373  include "pypars.inc"
1374 
1375  real forgetp,forintp
1376  external forgetp,forintp
1377  dimension sfm(8)
1378 
1379  t=y+rr*ta
1380  if(ita.eq.1 .or. ita.eq.4 .or. ita.eq.5)then
1381 
1382 * DIS case
1383  tt=sx-rr
1384  amf2=tt-t+amp2
1385  aks=t/tt
1386  anu=tt/ap
1387  epsi=ap2/tt
1388 
1389 * get unpolarized cross sections
1390  itara=int(rtara+0.1)
1391  itarz=int(rtarz+0.1)
1392 C....elke the order to call first mkr and than mkf2 is important for pythia6
1393 C to have R in mkf2 availabale
1394  call mkr(t,aks,dr)
1395 C if(aks.gt.0.99D0) then
1396 C write(*,*)"sx = ",sx , "y =", y, "ta = ", ta, "w2 =", w2
1397 C write(*,*)"q2 =",t, "xbj =",aks, "rr =", rr, "ita =", ita
1398 C endif
1399  call mkf2(t,aks,itara,itarz,f2,ff1)
1400  if (mstp(199).eq.0) then
1401  f1=amp*(1.+anu**2/t)*f2/anu/(1.+dr)
1402  else
1403  f1=ff1
1404  endif
1405 
1406 C call mkf2(t,aks,itara,itarz,f2,ff1)
1407 C call mkr(t,aks,dr)
1408 C f1=amp*(1.+anu**2/t)*f2/anu/(1.+dr)
1409 
1410 * get polarized cross sections
1411  if ((plrun.ne.0).and.(pnrun.ne.0)) then
1412  call mkasym(t,aks,itara,itarz,asym1,asym2)
1413 c ilyichev g1=f1*asym1*fdilut(t,aks,itara)
1414 
1415 cilyichev definition g1 and g2 via asym1, asym2 and f1.
1416  ga=sqrt(t)/anu
1417  g1=(asym1+ga*asym2)/(1+ga**2)*f1*fdilut(t,aks,itara)
1418  g2=(asym2-ga*asym1)/ga/(1+ga**2)*f1*fdilut(t,aks,itara)
1419  b1=0.d0
1420  b2=0.d0
1421  b3=0.d0
1422  b4=0.d0
1423  elseif((plrun.eq.0).and.(abs(pnrun).eq.2)) then
1424  call mkasym(t,aks,itara,itarz,asym1,asym2)
1425  b1=-3./2.*f1*asym1
1426  b2=2.d0*aks*b1
1427  b3=0.d0
1428  b4=0.d0
1429 C do we need
1430 C g1=0.d0
1431 C g2=0.d0
1432  else
1433  b1=0.d0
1434  b2=0.d0
1435  b3=0.d0
1436  b4=0.d0
1437  g1=0.d0
1438  g2=0.d0
1439  endif
1440  goto 10
1441  endif
1442 c
1443 c rtarn,rtarz,rtara are n,z,a of the nucleus
1444 c
1445  epsi=ap2/t
1446  rtarn=rtara-rtarz
1447 c
1448 c tf is t in fermi**(-2)
1449 c
1450  tf=t/chbar**2
1451 c
1452 c gep,gmp are electric and magnetic form factors of proton
1453 c s.i.bilenkaya et al pisma zhetf 19(1974)613
1454 c
1455  call ffpro(t,gep,gmp)
1456 * add neutron electric and magnetic form factors
1457 C-----fit of g.hohler et al.,nucl.phys. b114(1976)505.
1458 C FIT 8.2
1459 c
1460  f1s= 0.71/(0.783**2+t)-0.64/(1.02**2+t)-0.13/(1.80**2+t)
1461  f2s=-0.11/(0.783**2+t)+0.13/(1.02**2+t)-0.02/(1.80**2+t)
1462 c
1463  f1v= 0.05/(1.210**2+t)-0.52/(2.45**2+t)+0.28/(2.95**2+t)
1464  f2v=-1.99/(1.210**2+t)+0.20/(2.45**2+t)+0.19/(2.95**2+t)
1465 c
1466  f1rho=(0.955+0.090/(1+t/0.355)**2)/(1+t/0.536)/2.
1467  f2rho=(5.335+0.962/(1+t/0.268) )/(1+t/0.603)/2.
1468 c
1469  f1v=f1v+f1rho
1470  f2v=f2v+f2rho
1471 c
1472  f1p=f1s+f1v
1473  f1n=f1s-f1v
1474 c
1475  f2p=f2s+f2v
1476  f2n=f2s-f2v
1477 c
1478 * keep original proton form factors and use only neutron from this fit
1479 * gep=f1p-t*f2p/4./amp2
1480 * gmp=f1p+f2p
1481  gen=f1n-t*f2n/4./amp2
1482  gmn=f1n+f2n
1483 
1484  if(ita.eq.2)then
1485  tau=t/4./amt**2
1486  tau1=1.+tau
1487  if(ire.eq.0) then
1488  f1=2.*amp2*tau*gmn**2
1489  f2=4.*amp2*tau*(gen**2+tau*gmn**2)/tau1
1490  g1=amp2*tau*2.*gmn*(gen+tau*gmn)/tau1
1491  g2=2.*amp2*tau**2*gmn*(gen-gmn)/tau1
1492  b1=0d0
1493  b2=0d0
1494  b3=0d0
1495  b4=0d0
1496  elseif(ire.eq.1)then
1497  f1=2.*amp2*tau*gmp**2
1498  f2=4.*amp2*tau*(gep**2+tau*gmp**2)/tau1
1499  g1=amp2*tau*2.*gmp*(gep+tau*gmp)/tau1
1500  g2=2.*amp2*tau**2*gmp*(gep-gmp)/tau1
1501  b1=0d0
1502  b2=0d0
1503  b3=0d0
1504  b4=0d0
1505  elseif (ire.eq.2)then
1506  call ffdeu(t,fcdeu,fmdeu,fqdeu)
1507  fc=fcdeu*amt
1508  fm=fmdeu*amt
1509  fq=fqdeu*amt
1510  f1=4./3.*tau*tau1*fm**2
1511  f2=4./9.*tau*(8.*tau**2*fq**2+6.*tau*fm**2+9.*fc**2)
1512  g1=-1./3.*tau*fm*(2.*tau*fq+6.*fc+3.*tau*fm)
1513  g2=-1./3.*tau**2*fm*(2.*tau*fq+6*fc-3.*fm)
1514  b1=2.*tau**2*fm**2
1515  b2=4.*fm**2*tau**2+
1516  + 16./3.*tau**3/tau1*(tau*fq+3.*fc-3.*fm)*fq
1517  b3=-4./3.*(3.*tau+2.)*fm**2*tau**2+
1518  + 16./9.*tau**3/tau1*(tau*fq+3.*fc-3.*fm)*fq
1519  b4=4./3.*(6.*tau+1.)*fm**2*tau**2+
1520  + 16./9.*tau**3/tau1*(tau*fq+3.*fc+3.*(3.*tau+2.)*fm)*fq
1521  elseif (ire.eq.3)then
1522  call ffhe3(t,ge,gm)
1523  f1=2.*amt**2*tau*gm**2
1524  f2=4.*amt**2*tau*(ge**2+tau*gm**2)/tau1
1525  g1=amt**2*tau*2.*gm*(ge+tau*gm)/tau1
1526  g2=2.*amt**2*tau**2*gm*(ge-gm)/tau1
1527  b1=0d0
1528  b2=0d0
1529  b3=0d0
1530  b4=0d0
1531  elseif(ire.eq.4)then
1532  call ffhe4(t,ge,gm)
1533 *Test call ffhe3(t,ge,gm)
1534  f1=0d0
1535  f2=4.*amp2*tau*(rtarz*ge)**2
1536  g1=0d0
1537  g2=0d0
1538  b1=0d0
1539  b2=0d0
1540  b3=0d0
1541  b4=0d0
1542  elseif((ire.eq.14).or.(ire.eq.84).or.(ire.eq.20))then
1543  ff=forgetp(sngl(t))
1544 c if (t.lt.0.002) then
1545 c ff=forgetp(sngl(t))
1546 c print *,t,forgetp(sngl(t))
1547 c stop
1548 c endif
1549 * copy from Polrad - abr
1550  f1=0d0
1551  f2=4.*amp2*tau*(rtarz*ff)**2
1552  g1=0d0
1553  g2=0d0
1554  b1=0d0
1555  b2=0d0
1556  b3=0d0
1557  b4=0d0
1558  endif
1559  elseif (ita.eq.3) then
1560  tau=t/4./amp**2
1561  tau1=1.+tau
1562  call ffquas(t,geun,gmun,gepo,gmpo)
1563  f1=2.*amp2*tau*gmun**2
1564  f2=4.*amp2*tau*(geun**2+tau*gmun**2)/tau1
1565  g1=amp2*tau*2.*gmpo*(gepo+tau*gmpo)/tau1
1566  g2=2.*amp2*tau**2*gmpo*(gepo-gmpo)/tau1
1567  b1=0.
1568  b2=0.
1569  b3=0.
1570  b4=0
1571  endif
1572  10 continue
1573 
1574  sfm(1)=f1+qn/6.*b1
1575  sfm(2)=epsi*(f2+qn/6.*b2)
1576  sfm(3)=epsi*(g1+g2)
1577  sfm(4)=epsi**2*g2
1578  sfm(5)=epsi**2*b1
1579  sfm(6)=epsi**3*(b2/3.+b3+b4)
1580  sfm(7)=epsi*(b2/3.-b3)
1581  sfm(8)=epsi**2*(b2/3.-b4)
1582 
1583  return
1584  end
1585 
1586 C **********************************************************************
1587 C
1588 C FFPRO
1589 C
1590 C **********************************************************************
1591 
1592  subroutine ffpro(t,gep,gmp)
1593  implicit real*8(a-h,o-z)
1594 
1595  include "cmpcom.inc"
1596 
1597  gep=1.2742/(1.+t/0.6394**2)-.2742/(1.+t/1.582**2)
1598  gmp=(1.3262/(1.+t/0.6397**2)-.3262/(1.+t/1.3137**2))*amm
1599 c gep=1./((1.+.61*t)*(1.+2.31*t)*(1.+.04*t))
1600 c gmp=amm*gep
1601 
1602  end
1603 
1604 C **********************************************************************
1605 C
1606 C FFDEU
1607 C
1608 C **********************************************************************
1609 
1610  subroutine ffdeu(t,gc,gm,gq)
1611  implicit real*8(a-h,o-z)
1612 
1613  parameter(c2i3 = 2/3)
1614  dimension a(4),b(4),c(4)
1615  dimension al2ar(4),be2ar(4),ga2ar(4)
1616  data amd/1.8756280d0/chbar/0.19732858d0/
1617  data amp/0.9382796d0/
1618  data dmu/0.857406d0/dqu/25.84d0/
1619  integer para
1620 
1621  para=2
1622 
1623 C gd2=1d0/(1.d0+t/4./0.8952**2)**4
1624  gd2=1d0/(1.d0+t/4./0.89852**2)**4
1625  eta=t/4.d0/amd**2
1626  gd2e=gd2/(2d0*eta+1d0)
1627  sq2e=sqrt(2d0*eta)
1628 
1629  if(para.eq.1) then
1630  al2ar(1)=1.8591*chbar**2
1631  al2ar(4)=2d0*amd*0.58327d0
1632  be2ar(1)=19.586*chbar**2
1633  be2ar(4)=2d0*amd*0.1d0
1634  ga2ar(1)=1.0203*chbar**2
1635  ga2ar(4)=2d0*amd*0.17338d0
1636  else
1637  al2ar(1)=1.5250*chbar**2
1638  al2ar(4)=2d0*amd*0.24086d0
1639  be2ar(1)=43.677*chbar**2
1640  be2ar(4)=2d0*amd*0.029138d0
1641  ga2ar(1)=1.8706*chbar**2
1642  ga2ar(4)=2d0*amd*0.42693d0
1643  endif
1644 
1645  do i=2,3
1646  al2ar(i)=al2ar(1)+(al2ar(4)-al2ar(1))/3d0*(i-1)
1647  be2ar(i)=be2ar(1)+(be2ar(4)-be2ar(1))/3d0*(i-1)
1648  ga2ar(i)=ga2ar(1)+(ga2ar(4)-ga2ar(1))/3d0*(i-1)
1649  enddo
1650 
1651  if(para.eq.1) then
1652  a(1)=2.4818*chbar**2
1653  a(2)=-10.850*chbar**2
1654  a(3)=6.4416*chbar**2
1655  else
1656  a(1)=1.5706*chbar**2
1657  a(2)=12.238*chbar**2
1658  a(3)=-42.046*chbar**2
1659  endif
1660  a(4)=al2ar(4)*(1d0-a(2)/al2ar(2)-a(3)/al2ar(3)-a(1)/al2ar(1))
1661 
1662  if(para.eq.1) then
1663  b(1)=-1.7654*chbar
1664  b(2)=6.7874*chbar
1665  else
1666  b(1)=0.0704*chbar
1667  b(2)=0.1444*chbar
1668  endif
1669 
1670  bzn=1d0/be2ar(4)-1d0/be2ar(3)
1671  bbb=(2d0-dmu*amd/amp)/2d0/sqrt(2d0)/amd
1672  b(3)=(b(1)/be2ar(1)+b(2)/be2ar(2)-b(1)/be2ar(4)-b(2)/be2ar(4)
1673  + -bbb)/bzn
1674  b(4)=-(b(1)/be2ar(1)+b(2)/be2ar(2)-b(1)/be2ar(3)-b(2)/be2ar(3)
1675  + -bbb)/bzn
1676  ccc=(1d0-dmu*amd/amp-dqu)/4./amd**2
1677 
1678  if(para.eq.1) then
1679  c(1)=-0.053830d0
1680  else
1681  c(1)=-0.165770d0
1682  endif
1683 
1684  znc2=ga2ar(1)*(ga2ar(3)*ga2ar(4)-ga2ar(2)*ga2ar(3)
1685  + +ga2ar(2)**2-ga2ar(2)*ga2ar(4))
1686  c(2)=-ga2ar(2)/znc2*(c(1)*(
1687  + ga2ar(3)*ga2ar(4)-ga2ar(1)*ga2ar(3)+ga2ar(1)**2
1688  + -ga2ar(1)*ga2ar(4))-ccc*ga2ar(3)*ga2ar(4)*ga2ar(1) )
1689  znc3=ga2ar(1)*(ga2ar(3)-ga2ar(2))*(ga2ar(4)-ga2ar(3))
1690  c(3)=ga2ar(3)/znc3*(c(1)*(
1691  + ga2ar(2)*ga2ar(4)-ga2ar(1)*ga2ar(4)+ga2ar(1)**2
1692  + -ga2ar(1)*ga2ar(2))-ccc*ga2ar(2)*ga2ar(4)*ga2ar(1) )
1693  znc4=ga2ar(1)*(ga2ar(4)-ga2ar(2))*(ga2ar(4)-ga2ar(3))
1694  c(4)=-ga2ar(4)/znc4*(c(1)*(
1695  +ga2ar(2)*ga2ar(3)-ga2ar(1)*ga2ar(3)+ga2ar(1)**2-ga2ar(1)*ga2ar(2)
1696  + )-ccc*ga2ar(2)*ga2ar(3)*ga2ar(1) )
1697 
1698  g00=0d0
1699  gp0=0d0
1700  gpm=0d0
1701  sqt=sqrt(t)
1702  do i=1,4
1703  g00=g00+a(i)/(al2ar(i)+t)
1704  gp0=gp0+sqt*b(i)/(be2ar(i)+t)
1705  gpm=gpm+t*c(i)/(ga2ar(i)+t)
1706  enddo
1707 
1708  gc=gd2e*( (1d0-c2i3*eta)*g00+4d0*c2i3*sq2e*gp0
1709  + +c2i3*(2d0*eta-1d0)*gpm)
1710  gm=gd2e*(2d0*g00+2d0*(2d0*eta-1d0)/sq2e*gp0-2d0*gpm)
1711  gq=gd2e*(-g00+2d0/sq2e*gp0-(1d0+1d0/eta)*gpm)
1712 
1713  end
1714 
1715 C **********************************************************************
1716 C
1717 C FFHE3
1718 C
1719 C **********************************************************************
1720 
1721  subroutine ffhe3(t,ge,gm)
1722 c
1723 c l.i.shiff phys. rev. 133(1964)3b,802
1724 c
1725  implicit real*8(a-h,o-z)
1726 
1727  include "cmpcom.inc"
1728 
1729  tf=t/chbar**2
1730  qf=sqrt(tf)
1731  a=.675
1732  b=.366
1733  c=.836
1734  am=.654
1735  bm=.456
1736  cm=.821
1737  d=-6.78d-3
1738  p=.9
1739  q0=3.98
1740  f0=ddexp(-a**2*qf**2) - b**2*qf**2*ddexp(-c**2*qf**2)
1741  fm=ddexp(-am**2*qf**2) - bm**2*qf**2*ddexp(-cm**2*qf**2)
1742  df=d*ddexp(-((qf-q0)/p)**2)
1743  ge=f0+df
1744  gm=fm*rtara/rtarz * (-2.13)
1745 
1746  end
1747 C
1748  subroutine ffhe4(dq2,ge,gm)
1749 
1750 C charge form factor by SOG(Sum-of-Gaussian)
1751 C Helium-4,Carbon-12 : I.Sick, Phys.Lett.116B(1982)212.
1752 C For the Method see also I.Sick, Nucl.Phys.A218(1974)509.
1753 C or M.Arneodo, thesis Princeton 1991 page 270
1754 C Data: H. De Vries et al., Atomic Data and Nuclear Data Tables, 36 (1987) 495
1755 C Q2 IN 1/FERMI**2
1756 C DGAM2=(RP/SQRT(3/2))**2 = 0.66666666 for He4
1757  implicit real*8(a-h,o-z)
1758 
1759  include "cmpcom.inc"
1760 
1761  dimension dr(20),dq(20)
1762  DATA dr/.2d0,.6d0,.9d0,1.4d0,1.9d0,2.3d0,2.6d0
1763  : ,3.1d0,3.5d0,4.2d0,4.9d0,9*0.d0/
1764  DATA dq/.034724d0,.430761d0,.203166d0,.192986d0,.083866d0
1765  : ,.033007d0,.014201d0,0.d0,.00686d0,0.d0,.000438d0,9*0.d0/
1766  DATA dgam2/.666666666667d0/
1767  DATA nsog /11/
1768 
1769 * write(6,*) 'q2=',dq2
1770 
1771  dq2f=dq2/chbar**2
1772  dq1f=dsqrt(dq2f)
1773  delast=0.d0
1774  DO i=1,nsog
1775  d2rg=2.d0*dr(i)**2/dgam2
1776  dqri=dq1f*dr(i)
1777  dj0=1.d0
1778  IF(dqri.GE.1.d-6) dj0=dsin(dqri)/dqri
1779  dai=dq(i)/(1.d0+d2rg)*(dcos(dqri)+d2rg*dj0)
1780  delast=delast+dai
1781  enddo
1782  delast=dexp(-dq2f*dgam2/4.d0) * dabs(delast)
1783  ge=delast
1784  gm=0.0
1785 
1786  end
1787 
1788 
1789  REAL FUNCTION forgetp(Q2)
1790 C GET FORM FACTOR BY INTERPOLATION OUT OF TABLE FFNUC
1791 
1792  implicit none
1793  include "tailcom.inc"
1794  integer iq
1795  real q2
1796 
1797  forgetp=0.
1798 C LOWER BIN OF Q2
1799  iq=max0(int(q2/q2bin)+1,1)
1800  IF(iq+1.LT.nqbin) forgetp=ffnuc(iq)+
1801  1 (ffnuc(iq+1)-ffnuc(iq))*(q2/q2bin-float(iq-1))
1802 c if (iq.lt.20) print *,iq,ffnuc(iq),forgetp
1803 c print *,iq,q2,ffnuc(iq),forgetp
1804  RETURN
1805  END
1806 
1807  SUBROUTINE fordop
1808 C produce form factor table by fourier-transformation
1809 C of the charge distribution
1810 
1811  implicit none
1812  include "tailcom.inc"
1813  real q2max,rmin,rmax,eps,hbc
1814  real q2
1815  real forintp,gauss
1816  external forintp
1817  integer nqbin1,iq
1818 
1819  DATA q2max,rmin,rmax,eps,nqbin1/.3,1.e-6,20.,1.e-6,600/
1820  DATA hbc /197.3234/
1821 C
1822  nqbin=nqbin1
1823 C
1824 C loop over q-bins
1825  q2bin=q2max/float(nqbin-1)
1826 C
1827  DO 10 iq=1,nqbin
1828 C first bin is q2=0. last is q2=q2max
1829  q2=q2bin*float(iq-1)
1830 C calc 3-momentum forq (transform from gev to fermi?)
1831  qfor=sqrt(q2)*1000./hbc
1832 C perform fourier-transform by integ with CERN-Gauss routine
1833  ffnuc(iq)=gauss(forintp,rmin,rmax,eps)
1834 C normalise form factor in a way, that FFNUC(Q=0)=1.
1835  IF(iq.GT.1)ffnuc(iq)=ffnuc(iq)/ffnuc(1)
1836 c if (iq.lt.20) print *,iq,q2,ffnuc(iq)
1837 10 CONTINUE
1838  ffnuc(1)=1.
1839 C
1840  RETURN
1841  END
1842 
1843 
1844  REAL FUNCTION forintp(R)
1845 C INTEGRAND OF FOURIER TRANSFORMATION INTEGRAL
1846 C THREE-DIMENSIONAL INTEGRATION IS REDUCED TO
1847 C AN INTEGRATION OVER THE RADIUS BY MULTIPLYING
1848 C WITH A BESSEL FKT.
1849 
1850  implicit none
1851  include "tailcom.inc"
1852  real r,qr,forrhop
1853 
1854  IF (qfor.GT.1.e-5) THEN
1855  qr=qfor*r
1856  forintp=forrhop(r)*sin(qr)*r**2/qr
1857  ELSE
1858  forintp=forrhop(r)*r**2
1859  ENDIF
1860 c print *,'forintp=',forintp
1861  RETURN
1862  END
1863 
1864 
1865  REAL FUNCTION forrhop(R)
1866 C charge distribution as fkt. of radius
1867 C generalized 2-parameter fermi for not too light nuclei
1868 C Ref: T.de Forest and J.D.Walecka, Adv.in Phys.15(1966)1 -- page 21
1869 C they cite R.Hofstadter, Rev Mod Phys 28 (1963) 214
1870 C z = 2.4 / (4*ln(3))
1871 
1872  implicit none
1873  include "tailcom.inc"
1874  include "cmpcom.inc"
1875  real r,zpara,cpara,wpara
1876 
1877  if (ire.eq.14) then
1878 C - 3 parameter Fermi charge distr of 14N
1879 C H. de Vries et al, Atomic Data and Nuclear Tables 36, 495-536 (1987)
1880  cpara=2.57
1881  zpara=0.5052
1882  wpara=-0.180
1883  forrhop=(1.0+wpara*(r/cpara)**2)/
1884  1 (1.0+exp((r-cpara)/zpara))
1885  elseif (ire.eq.20) then
1886 C - 3 parameter Fermi charge distr of 20 Ne
1887 C H. de Vries et al, Atomic Data and Nuclear Tables 36, 495-536 (1987)
1888  cpara=2.791
1889  zpara=0.698
1890  wpara=-0.168
1891  forrhop=(1.0+wpara*(r/cpara)**2)/
1892  1 (1.0+exp((r-cpara)/zpara))
1893  elseif (ire.eq.84) then
1894 C - RHO2BOMO (TERADNMC.CAR)
1895 C generalized 2-parameter fermi
1896 C Ref: Bohr & Mottelson, Nuclear Stucture (Benjamin, New York 1969)
1897 C cited in Date et al., PRL 52 (1984) 2344
1898  zpara=0.54
1899  cpara=1.12*rtara**(1./3.)-0.86/rtara**(1./3.)
1900  forrhop=1./(1.+exp((r-cpara)/zpara))
1901  else
1902 C generalized 2-parameter fermi for not too light nuclei
1903 C Ref: T.de Forest and J.D.Walecka, Adv.in Phys.15(1966)1 -- page 21
1904 C they cite R.Hofstadter, Rev Mod Phys 28 (1963) 214
1905 C z = 2.4 / (4*ln(3))
1906  zpara=0.546144
1907  cpara=1.07*rtara**(1./3.)
1908  forrhop=1./(1.+exp((r-cpara)/zpara))
1909  endif
1910  RETURN
1911  END
1912 
1913 
1914 C **********************************************************************
1915 C
1916 C FFQUAS
1917 C
1918 C **********************************************************************
1919 
1920  subroutine ffquas(t,geun,gmun,gepo,gmpo)
1921  implicit real*8(a-h,o-z)
1922 
1923  include "cmpcom.inc"
1924  include "tailcom.inc"
1925 
1926  call ffpro(t,gep,gmp)
1927  tf=t/chbar**2
1928  tau=t/4./amp**2
1929  tau1=1.+tau
1930 c
1931 c t. de forest and d.j. valecka adv. in phys. 15(1966) no.57
1932 c
1933  if(ire.eq.3)then
1934  supele=1.
1935  supmag=1.
1936  if(t.le.(2.d0*fermom)**2)then
1937  sqrat=dsqrt(t)/fermom
1938  supele=0.75*sqrat-sqrat**3/16.
1939  supmag=supele
1940  endif
1941  else
1942  supele=supst(t)
1943 c qbold=sqrt(t*tau1)
1944 c supele=1.-dsbern(qbold)**2
1945  supmag=supele
1946  endif
1947  geun=gep*dsqrt(supele*rtarz)
1948  rtarn=rtara-rtarz
1949  gmun=gep*dsqrt(supmag*(rtarz*amm**2+rtarn*amn**2))
1950  if(ire.eq.2)then
1951  gepo=geun
1952  gmpo=gmun
1953  else
1954  gepo=0.
1955  rtarn=rtara-rtarz
1956  gmpo=gep*dsqrt(supmag*(rtarn*amn**2))
1957  endif
1958  end
1959 
1960 ********************** supst ************************************
1961 
1962  double precision function supst(t)
1963  implicit real*8(a-h,o-z)
1964  data chbar/.197328d0/
1965 c
1966 c tf is t in fermi**(-2)
1967 c
1968  tf=t/chbar**2
1969 c
1970 c s.stein et al phys. rev. 12(1975)1884 (appendix 1)
1971 c
1972  sqtf=dsqrt(tf)
1973  delff=(datan(sqtf/.93d0)-2.*datan(sqtf/3.19d0)+
1974  + datan(sqtf/5.45d0))*1.580d0/sqtf
1975  delff=dmax1(0.d0,delff)
1976  supele=1.-delff**2
1977  supst=dmax1(0.d0,supele)
1978  return
1979  end
1980 
1981 C **********************************************************************
1982 C
1983 C DDEXP
1984 C
1985 C **********************************************************************
1986 
1987  double precision function ddexp(x)
1988  implicit real*8(a-h,o-z)
1989 
1990  ddexp=0.
1991  if(x.gt.-50.)ddexp=exp(x)
1992 
1993  return
1994  end
1995 
1996 C **********************************************************************
1997 C
1998 C DQN32
1999 C
2000 C **********************************************************************
2001 
2002  subroutine dqn32(xl,xu,fct,y)
2003 c
2004 c
2005  external fct
2006  double precision xl,xu,y,a,b,c,fct
2007 c
2008  a=.5d0*(xu+xl)
2009  b=xu-xl
2010  c=.49863193092474078d0*b
2011  y=.35093050047350483d-2*(fct(a+c)+fct(a-c))
2012  c=.49280575577263417d0*b
2013  y=y+.8137197365452835d-2*(fct(a+c)+fct(a-c))
2014  c=.48238112779375322d0*b
2015  y=y+.12696032654631030d-1*(fct(a+c)+fct(a-c))
2016  c=.46745303796886984d0*b
2017  y=y+.17136931456510717d-1*(fct(a+c)+fct(a-c))
2018  c=.44816057788302606d0*b
2019  y=y+.21417949011113340d-1*(fct(a+c)+fct(a-c))
2020  c=.42468380686628499d0*b
2021  y=y+.25499029631188088d-1*(fct(a+c)+fct(a-c))
2022  c=.39724189798397120d0*b
2023  y=y+.29342046739267774d-1*(fct(a+c)+fct(a-c))
2024  c=.36609105937014484d0*b
2025  y=y+.32911111388180923d-1*(fct(a+c)+fct(a-c))
2026  c=.33152213346510760d0*b
2027  y=y+.36172897054424253d-1*(fct(a+c)+fct(a-c))
2028  c=.29385787862038116d0*b
2029  y=y+.39096947893535153d-1*(fct(a+c)+fct(a-c))
2030  c=.25344995446611470d0*b
2031  y=y+.41655962113473378d-1*(fct(a+c)+fct(a-c))
2032  c=.21067563806531767d0*b
2033  y=y+.43826046502201906d-1*(fct(a+c)+fct(a-c))
2034  c=.16593430114106382d0*b
2035  y=y+.45586939347881942d-1*(fct(a+c)+fct(a-c))
2036  c=.11964368112606854d0*b
2037  y=y+.46922199540402283d-1*(fct(a+c)+fct(a-c))
2038  c=.7223598079139825d-1*b
2039  y=y+.47819360039637430d-1*(fct(a+c)+fct(a-c))
2040  c=.24153832843869158d-1*b
2041  y=b*(y+.48270044257363900d-1*(fct(a+c)+fct(a-c)))
2042 
2043  return
2044  end
2045 
2046 C **********************************************************************
2047 C
2048 C INTBTK2
2049 C
2050 C **********************************************************************
2051 
2052  subroutine intbtk2(dbtk,nbtk,dtkmax)
2053  implicit real*8 (a-h,o-z)
2054 c
2055 c *** choose 7 bins for integration over the theta k angle
2056  dimension dbtk(8)
2057 
2058  include "sxycom.inc"
2059  include "cmpcom.inc"
2060  include "ppicom.inc"
2061  include "radgen.inc"
2062 
2063  data dc /3.5d0/
2064 c
2065 c bins are chosen according to the s peak and p peak
2066 c width of the peaks:
2067 
2068  dwids=dsqrt(ap*aml/s)/dc
2069  dwidp=dsqrt(ap*aml/x)/dc
2070 * avoid values above 1.0 (for x values below 1e-07, abr 19.02.04)
2071 * dts=acos( (s*sx+ap2*y)/sqls/sqly )
2072 * dtp=acos( (x*sx-ap2*y)/sqlx/sqly )
2073  dts=acos( min(1.0d0,(s*sx+ap2*y)/sqls/sqly) )
2074  dtp=acos( min(1.0d0,(x*sx-ap2*y)/sqlx/sqly) )
2075 
2076  dpi=pi
2077 
2078 c
2079 c define bin boundaries
2080  dbtk(1)=0.d0
2081  dbtk(2)=dmin1(4.d0*dc*dwids,dts/3.0d0)
2082  dbtk(3)=dmax1(dts-dwids,dts/1.5d0)
2083  dbtk(4)=dmin1(dts+dwids,dts+(dtp-dts)/3.0d0)
2084  dbtk(5)=dmax1(dtp-dwidp,dts+(dtp-dts)/1.5d0)
2085  dbtk(6)=dmin1(dtp+dwidp,dtp+(dpi-dtp)/3.0d0)
2086  dbtk(7)=.98d0*dbtk(6)+.02d0*dpi
2087  dbtk(8)=dpi
2088 c
2089 c check cut in theta which follows from infra red cut:
2090 c
2091  nbtk=0
2092  do 10 ka=2,8
2093  nbtk=nbtk+1
2094  if(dbtk(nbtk+1).ge.dtkmax) goto 20
2095  10 continue
2096  goto 30
2097  20 dbtk(nbtk+1)=dmin1(dbtk(nbtk+1),dtkmax)
2098  30 continue
2099 c
2100 c nbtk=number of valid bins
2101 c
2102  return
2103  end
2104 
2105 C **********************************************************************
2106 C
2107 C INTTK2
2108 C
2109 C **********************************************************************
2110 
2111  subroutine inttk2(isumtk,dbtk1,dbtk2,dsumtk,derrtk)
2112  implicit real*8 (a-h,o-z)
2113 
2114 c integrate over ntk bins of brems-gamma theta angle dtk
2115  include "radgen.inc"
2116  include "sxycom.inc"
2117  include "cmpcom.inc"
2118  include "tailcom.inc"
2119  include "ppicom.inc"
2120 
2121  ddtk=(dbtk2-dbtk1)/ntk
2122  dtk=dbtk1-ddtk*.5d0
2123  dsum=0.d0
2124  derr1=0.d0
2125 c loop over all bins
2126  do 10 itk=1,ntk
2127  isumtk=isumtk+1
2128  dtk=dtk+ddtk
2129  itkcur=isumtk
2130 c calculate integrand
2131 C call intsai(dtk,dtsai)
2132 C dsum=dsum+dtsai
2133  ta=(sx-sqly*cos(dtk))/ap2
2134  dsigma=an*alfa/pi*rv2(ta) * sin(dtk)*sqly/ap2
2135  dsum=dsum+dsigma
2136 C write(9,*)dsigma,ta,dtk
2137  if(ita.le.3)then
2138  dsitkm(isumtk,ita) = dsumtk + ddtk * dsum
2139  dcvtkm(isumtk,ita) = dtk
2140  ddetkm(isumtk,ita) = ddtk
2141  endif
2142 c add up errors:
2143 c errors of theta integration:
2144  if (itk.gt.2.and.itk.lt.ntk) then
2145  derr1=derr1+dabs(dsigma-dold)
2146  elseif(itk.eq.2.or.itk.eq.ntk) then
2147  derr1=derr1+1.5d0*dabs(dsigma-dold)
2148  endif
2149 c
2150  dold=dsigma
2151 c
2152  10 continue
2153 c
2154 c integral on big bin is:
2155 C print *,dsumtk,ddtk,dsum
2156  dsumtk=dsumtk+ddtk*dsum
2157  derrtk=derrtk+ddtk*derr1
2158 c
2159  return
2160  end
2161 
2162 
2163 C **********************************************************************
2164 C
2165 C ITAFUN
2166 C
2167 C **********************************************************************
2168 
2169  subroutine itafun(ys,xs,pl,pn,ixytest,itagen)
2170 
2171  include "cmpcom.inc"
2172  include "radgen.inc"
2173  include "density.inc"
2174  include "xytabcom.inc"
2175 
2176  dimension xym(2),na(2),a(2*nt)
2177 
2178  do nny=nty,1,-1
2179  if(ys.gt.y(nny))goto 10
2180  enddo
2181  5 print *,' ys=',ys,' out of the region'
2182  print *, 'y(nny)=', y(nny), nny, nt, nty
2183  stop
2184  10 if(nny.eq.nt)goto 5
2185 
2186  if(ixytest.eq.1)then
2187  stop .eq.' ixytest1 is not supported in the version'
2188  else
2189  xym(1)=xs
2190  xym(2)=ys
2191  na(1)=ntx
2192  na(2)=nty
2193  do i=1,ntx
2194  a(i)=x(i)
2195  a(i+ntx)=y(i)
2196  enddo
2197 * what nonsense if this ? - abr 19.02.04
2198 * none of the polarised quantities (tbor_p etc) is ever different from zero !
2199 * use final quantities (tbor etc) immediately
2200 * different array sizes for DIS and phoitoproduction kinematics
2201 * also include other quantities as redfac to be available with LUT
2202  if (ntx.eq.ntdis) then
2203  tbor=fint(2,xym,na,a,tbor_udis)
2204  sigrad=fint(2,xym,na,a,sigrad_udis)
2205  sig1g=fint(2,xym,na,a,sig1g_udis)
2206  vac=fint(2,xym,na,a,vac_udis)
2207  vertex=fint(2,xym,na,a,vertex_udis)
2208  small=fint(2,xym,na,a,small_udis)
2209  redfac=fint(2,xym,na,a,redfac_udis)
2210  else
2211  tbor=fint(2,xym,na,a,tbor_upho)
2212  sigrad=fint(2,xym,na,a,sigrad_upho)
2213  sig1g=fint(2,xym,na,a,sig1g_upho)
2214  vac=fint(2,xym,na,a,vac_upho)
2215  vertex=fint(2,xym,na,a,vertex_upho)
2216  small=fint(2,xym,na,a,small_upho)
2217  redfac=fint(2,xym,na,a,redfac_upho)
2218  endif
2219 * tbor=tboru+pl*pn*tborp
2220 * sigrad=sigradu+pl*pn*sigradp
2221 * sig1g=sig1gu+pl*pn*sig1gp
2222 
2223  endif
2224 
2225  sigcor=sigrad/sig1g
2226 !|bs>
2227 cc write(*,*)'radgen:',xs,ys,sigcor,sigrad,sig1g,sig1gu,pl,pn,sig1gp
2228 
2229  r1=rlu(0)
2230  rr1=r1*sigrad
2231 
2232 * same nonsense here (no need for tineu and tinep, use tine immediately)
2233  if (ntx.eq.ntdis) then
2234  tine=fint(2,xym,na,a,tine_udis)
2235  tnuc=fint(2,xym,na,a,tnuc_udis)
2236  tpro=fint(2,xym,na,a,tpro_udis)
2237  else
2238  tine=fint(2,xym,na,a,tine_upho)
2239  tnuc=fint(2,xym,na,a,tnuc_upho)
2240  tpro=fint(2,xym,na,a,tpro_upho)
2241  endif
2242 
2243  if(rr1.gt.sigrad-tbor)then
2244  itagen=0
2245  return
2246  endif
2247 
2248 C write(*,'(6g13.5)')sigrad,tbor,tine,tpro,tnuc,rr1
2249 c pause
2250 
2251  if(rr1.gt.(tpro+tnuc)) then
2252  itagen=1
2253 c scgen=rr1-tpro-tnuc
2254 
2255  elseif(rr1.gt.tnuc)then
2256  itagen=3
2257 c scgen=rr1-tnuc
2258  else
2259 c scgen=rr1
2260  itagen=2
2261  endif
2262 
2263 ! Speed optimization by Joe Seele (seele@uiuc.edu) June 4, 2002
2264 
2265 C #undef OLDCODE
2266  do ii=1,7*ntk
2267 C #ifdef OLDCODE
2268 C do ix=1,ntx
2269 C do iy=1,nty
2270 C if (ntx.le.ntdis) then
2271 C wrk(ix,iy)=densdis(ntx+1-ix,iy,ii,itagen)
2272 C else
2273 C wrk(ix,iy)=denspho(ntx+1-ix,iy,ii,itagen)
2274 C endif
2275 C enddo
2276 C enddo
2277 C dsitkm(ii,itagen)=fint(2,xym,na,a,wrk)
2278 C #else
2279  if (ntx.le.ntdis) then
2280  dsitkm(ii,itagen)=fint(2,xym,na,a,densdis(1,1,ii,itagen))
2281  else
2282  dsitkm(ii,itagen)=fint(2,xym,na,a,denspho(1,1,ii,itagen))
2283  endif
2284 C #endif
2285  enddo
2286 
2287  do ii=1,7
2288 C #ifdef OLDCODE
2289 C do ix=1,ntx
2290 C do iy=1,nty
2291 C if (ntx.le.ntdis) then
2292 C wrk(ix,iy)=widdis(ntx+1-ix,iy,ii,itagen)
2293 C else
2294 C wrk(ix,iy)=widpho(ntx+1-ix,iy,ii,itagen)
2295 C endif
2296 C enddo
2297 C enddo
2298 C ddetkm(ii,itagen)=fint(2,xym,na,a,wrk)
2299 C #else
2300  if (ntx.le.ntdis) then
2301  ddetkm(ii,itagen)=fint(2,xym,na,a,widdis(1,1,ii,itagen))
2302  else
2303  ddetkm(ii,itagen)=fint(2,xym,na,a,widpho(1,1,ii,itagen))
2304  endif
2305 C #endif
2306  enddo
2307 
2308  do iittkk=1,7*ntk
2309  ssuumm=0.
2310  do itoo=1,iittkk/ntk
2311  ssuumm=ssuumm+ddetkm(itoo,itagen)*ntk
2312  enddo
2313  dcvtkm(iittkk,itagen)=ssuumm+ddetkm((iittkk-1)/ntk+1,itagen)
2314  + *(mod(iittkk,ntk)-0.5)
2315  enddo
2316 
2317  end
2318 
2319 
2320 
2321 C **********************************************************************
2322 C
2323 C XYTABL
2324 C
2325 C **********************************************************************
2326 
2327  subroutine xytabl(tname,e1,plrun,pnrun,ixytest,ire)
2328 
2329  include "cmpcom.inc"
2330  include "radgen.inc"
2331  include "density.inc"
2332  include "xytabcom.inc"
2333  include "mc_set.inc"
2334 
2335  character*256 tname
2336 
2337  data lun/66/
2338  open (unit=lun,file=tname)
2339 
2340  if(ixytest.eq.0)then
2341 
2342  do iy=1,nty
2343  y(iy)=radgen_ymin+(radgen_ymax-radgen_ymin)/(nty-1)*(iy-1)
2344  enddo
2345 
2346 *abr step=(radgen_xmax-radgen_xmin)/(nt-1)
2347 * abr 19.02.04
2348  do ix=1,ntx
2349 * start from high x, then bring in ascending order
2350  if (ix.le.18) then
2351  step=0.05
2352  x(ix)=1.0-step*ix
2353  elseif (ix.le.27) then
2354  step=0.01
2355  x(ix)=x(18)-step*(ix-18)
2356  elseif (ix.le.36) then
2357  step=0.001
2358  x(ix)=x(27)-step*(ix-27)
2359  elseif (ix.le.48) then
2360  if (mod(ix,2).gt.0) then
2361  x(ix)=x(ix-1)/2
2362  else
2363  x(ix)=x(ix-1)/5
2364  endif
2365  endif
2366  enddo
2367  CALL flpsor(x,ntx)
2368 
2369  do iy=1,nty
2370  do ix=1,ntx
2371 *abr x(ix,iy)=radgen_xmin+step*(ix-1)
2372  bmp=amp
2373  bmc2=1.16
2374  ylim=(bmc2-bmp**2)/(2.*bmp*e1*(1d0-x(ix)))
2375  ys=y(iy)
2376  if(ys.lt.ylim ) ys=ylim
2377  call mpolrad(e1,ys,x(ix),1.,plrun,pnrun,-1)
2378  tbor_u(ix,iy)=tbor
2379  sig1g_u(ix,iy)=sig1g
2380  sigrad_u(ix,iy)=sigrad
2381  tine_u(ix,iy)=tine
2382  tnuc_u(ix,iy)=tnuc
2383  tpro_u(ix,iy)=tpro
2384  vac_u(ix,iy)=vac
2385  vertex_u(ix,iy)=vertex
2386  small_u(ix,iy)=small
2387  redfac_u(ix,iy)=redfac
2388 C #ifdef TEST
2389 C write(*,'(a3,6g12.4)')' u ',
2390 C + y(iy),x(ix),sig1g,tnuc,tpro,sigrad
2391 C #endif
2392  write(lun,'(1x,14g14.4)')x(ix),y(iy)
2393  + ,sig1g_u(ix,iy),sigrad_u(ix,iy),tbor_u(ix,iy)
2394  + ,tine_u(ix,iy),tnuc_u(ix,iy),tpro_u(ix,iy)
2395  + ,vac_u(ix,iy),vertex_u(ix,iy),small_u(ix,iy)
2396  + ,redfac_u(ix,iy)
2397 
2398  do itkm=0,6
2399  width(ix,iy,itkm+1,1)=ddetkm(itkm*ntk+1,1)
2400  width(ix,iy,itkm+1,2)=ddetkm(itkm*ntk+1,2)
2401  if(ire.ne.1)
2402  + width(ix,iy,itkm+1,3)=ddetkm(itkm*ntk+1,3)
2403  enddo
2404 
2405  ittmax=3
2406  if(ire.eq.1)ittmax=2
2407  write(lun,*)((width(ix,iy,itkm+1,itt)
2408  + ,itkm=0,6),itt=1,ittmax),ndxtkm
2409 
2410  do itkm=1,ntk*7
2411 * avoid values smaller than E-24, happens in case of N14 (abr)
2412 * similar problem seems to there at high x values
2413 * is this really the solution ???
2414 * at least it seems to work....
2415 * denstk(ix,iy,itkm,1)=dsitkm(itkm,1)/dsitkm(ntk*7,1)
2416  if (dsitkm(itkm,1).lt.1.e-24 .or.
2417  & dsitkm(ntk*7,1).lt.1e-24 ) then
2418  denstk(ix,iy,itkm,1)=1.0
2419  else
2420  denstk(ix,iy,itkm,1)=dsitkm(itkm,1)/dsitkm(ntk*7,1)
2421  endif
2422 
2423  if (dsitkm(itkm,2).lt.1.e-24 .or.
2424  & dsitkm(ntk*7,2).lt.1e-24 ) then
2425  denstk(ix,iy,itkm,2)=1.0
2426  else
2427  denstk(ix,iy,itkm,2)=dsitkm(itkm,2)/dsitkm(ntk*7,2)
2428  endif
2429 
2430 c if (ix.gt.20) then
2431 c write(6,*) dsitkm(itkm,2),dsitkm(ntk*7,2),
2432 c + denstk(ix,iy,itkm,2)
2433 c endif
2434 * 28.01.04 - better include safety limits here as well
2435 * if(ire.ne.1)
2436 * + denstk(ix,iy,itkm,3)=dsitkm(itkm,3)/dsitkm(ntk*7,3)
2437  if(ire.ne.1) then
2438  if (dsitkm(itkm,3).lt.1.e-24 .or.
2439  & dsitkm(ntk*7,3).lt.1e-24 ) then
2440  denstk(ix,iy,itkm,3)=1.0
2441  else
2442  denstk(ix,iy,itkm,3)=dsitkm(itkm,3)/dsitkm(ntk*7,3)
2443  endif
2444  endif
2445  enddo
2446  write(lun,'(35g14.4)')(denstk(ix,iy,itkm,1)
2447  + ,itkm=1,7*ntk)
2448  write(lun,'(35g14.4)')(denstk(ix,iy,itkm,2)
2449  + ,itkm=1,7*ntk)
2450  if(ire.ne.1)
2451  + write(lun,'(35g14.4)')(denstk(ix,iy,itkm,3)
2452  + ,itkm=1,7*ntk)
2453  enddo
2454  enddo
2455  close(lun)
2456 
2457  elseif(ixytest.gt.0)then
2458  do iy=1,nty
2459  do ix=1,ntx
2460 C #ifdef TEST
2461 C print *,iy,ix
2462 C #endif
2463  if (ntx.eq.ntdis) then ! DIS range (xmin=0.002)
2464  read(lun,*)x(ix),y(iy)
2465  + ,sig1g_udis(ix,iy),sigrad_udis(ix,iy),tbor_udis(ix,iy)
2466  + ,tine_udis(ix,iy),tnuc_udis(ix,iy),tpro_udis(ix,iy)
2467  + ,vac_udis(ix,iy),vertex_udis(ix,iy)
2468  + ,small_udis(ix,iy),redfac_udis(ix,iy)
2469 c print *,ix,iy
2470 c + ,sig1g_u(ix,iy),sigrad_u(ix,iy),tbor_u(ix,iy)
2471 c + ,tine_u(ix,iy),tnuc_u(ix,iy),tpro_u(ix,iy)
2472 
2473  ittmax=3
2474  if(ire.eq.1)ittmax=2
2475  read(lun,*)
2476  + ((widdis(ix,iy,ite,itt),ite=1,7),itt=1,ittmax),ndxtkm
2477  read(lun,*)(densdis(ix,iy,itkm,1)
2478  + ,itkm=1,7*ntk)
2479  read(lun,*)(densdis(ix,iy,itkm,2)
2480  + ,itkm=1,7*ntk)
2481  if(ire.ne.1)
2482  + read(lun,*)(densdis(ix,iy,itkm,3)
2483  + ,itkm=1,7*ntk)
2484  else ! photoproduction (xmin=1e-07)
2485  read(lun,*)x(ix),y(iy)
2486  + ,sig1g_upho(ix,iy),sigrad_upho(ix,iy),tbor_upho(ix,iy)
2487  + ,tine_upho(ix,iy),tnuc_upho(ix,iy),tpro_upho(ix,iy)
2488  + ,vac_upho(ix,iy),vertex_upho(ix,iy)
2489  + ,small_upho(ix,iy),redfac_upho(ix,iy)
2490 
2491  ittmax=3
2492  if(ire.eq.1)ittmax=2
2493  read(lun,*)
2494  + ((widpho(ix,iy,ite,itt),ite=1,7),itt=1,ittmax),ndxtkm
2495  read(lun,*)(denspho(ix,iy,itkm,1)
2496  + ,itkm=1,7*ntk)
2497  read(lun,*)(denspho(ix,iy,itkm,2)
2498  + ,itkm=1,7*ntk)
2499  if(ire.ne.1)
2500  + read(lun,*)(denspho(ix,iy,itkm,3)
2501  + ,itkm=1,7*ntk)
2502  endif
2503 
2504  enddo
2505  enddo
2506  else
2507  stop 'xytabl -- ixytest < 0'
2508  endif
2509  close(lun)
2510 
2511  end
2512 
2513 C **********************************************************************
2514 C
2515 C
2516 C
2517 C **********************************************************************