Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:21:21 

0001 
0002 C **********************************************************************
0003 C
0004 C     RADATA
0005 C
0006 C **********************************************************************
0007 
0008       BLOCK DATA RADATA
0009       implicit real*8(a-h,o-z)
0010 
0011       include "cmpcom.inc"
0012       include "ppicom.inc"
0013 
0014       data
0015      +     amm/2.7928456d0/,
0016      +     amn/-1.913148d0/,
0017      +     chbar/.197328d0/,
0018      +     barn/.389379d6/,
0019      +     aml/.511000d-3/,
0020      +     aml2/.261112d-6/,
0021      +     al2/.522240d-6/,
0022      +     isf20/8/,
0023      +     pi/3.1415926d0/,
0024      +     pi2/9.869604d0/,
0025      +     alfa/.729735d-2/,
0026      +     amc2/1.151857d0/,
0027      +     amp/.938272d0/,
0028      +     amh/.938272d0/,
0029      +     i2/1,1,1,2,3,3,1,2/,
0030      +     i1/3,3,4,4,3,3,3,3/
0031 
0032       end
0033 
0034 C **********************************************************************
0035 C
0036 C     RADGEN
0037 C
0038 C **********************************************************************
0039 
0040       SUBROUTINE RADGEN(e1,q2l,nul,ys,xs,phil,PhRAD,q2tr,anutr,WEIGHT)
0041       IMPLICIT NONE
0042 ***********************************************************
0043 *  INPUT :
0044 *        e1    : beam energy
0045 *        q2,nu : uncorrected kinematics of scattered lepton
0046 *        phil  : azimuthal angle of scattered lepton
0047 *  OUTPUT :
0048 *          PhRAD(4) : real rad. photon 4 vector
0049 *          q2tr     : true Q2
0050 *          Utr      : true U
0051 ***********************************************************
0052 
0053       include "cmpcom.inc"
0054       include "radgen.inc"
0055       include "phiout.inc"
0056       include "tailcom.inc"
0057       include "radgenkeys.inc"
0058 
0059       REAL PhRAD(*),q2tr,anutr
0060       REAL e1,q2l,nul,phil,xs,ys,weight
0061       INTEGER itagen
0062 *
0063 
0064 * reset proton mass (abr)
0065 * necessary to get correct kinematics after first elastic event !
0066       amp=0.938272d0
0067 
0068 C... Don't recalculate x and y, but use the ones already calculated
0069 C      ys=geny
0070 C      xs=genx
0071 C      ys=nul/e1
0072 C      xs=q2l/(2.*amp*nul)
0073 
0074       write(*,*)"radgen: ",ys, q2l, xs, nul, phil, e1
0075 C Do we want to use a lookup table?
0076       if(ixytest.ge.0)then
0077         if(kill_elas_res.ne.2)then
0078           call itafun(ys,xs,plrun,pnrun,ixytest,itagen)
0079         else
0080           itagen=0
0081           sigcor=1.
0082         endif
0083       else
0084         itagen=-1
0085       endif
0086 
0087       if(itagen.ne.0.and.abs(ixytest).ne.2)then
0088         iphi=0
0089         call mpolrad(e1,ys,xs,1.,plrun,pnrun,itagen)
0090         if(ixytest.eq.1)ixytest=-2
0091       endif
0092       call radgam_pol(e1,ys,xs,phil,ixytest,itagen,q2tr,anutr)
0093       
0094 * set kinematics of possibly radiated real gamma
0095       phrad(1)=sngl(dplabg(1))
0096       phrad(2)=sngl(dplabg(2))
0097       phrad(3)=sngl(dplabg(3))
0098 C to avoid rounding errors !
0099 CCC      phrad(4)=ys*e1-anutr
0100       phrad(4)=sqrt(phrad(1)*phrad(1)+phrad(2)*phrad(2)+
0101      +              phrad(3)*phrad(3))
0102 
0103 * set the weight
0104       if(kill_elas_res.eq.2)then
0105         weight=1.
0106       else
0107         weight = sigcor
0108       endif
0109 
0110       return
0111       end
0112 
0113 C **********************************************************************
0114 C
0115 C     PHIDST_POL
0116 C
0117 C **********************************************************************
0118 
0119       subroutine phidst_pol
0120 *     -----------------
0121 *
0122 *     calculate phi-distribution of bremsstrahlung photons
0123 *     needed for correction of hadronic distributions
0124 *
0125       include "double.inc"
0126       include "phiout.inc"
0127       include "cmpcom.inc"
0128       include "radgen.inc"
0129       include "ppicom.inc"
0130 C
0131       dimension db(3),dx(31),dy(31),dz(31)
0132       external dfufi_pol
0133       data ndim /31/
0134 C
0135 C     boundaries for 2 integrations
0136 C
0137       db(1)=0.d0
0138       db(2)=pi/10.d0
0139       db(3)=pi
0140       k=0
0141       dsum=0.d0
0142       do i=1,2
0143         dbnd1=db(i)
0144         dbnd2=db(i+1)
0145         call radquad_pol(dfufi_pol,dbnd1,dbnd2,dx,dy,dz,dh,ndim)
0146         do j=1,ndim
0147              if(i.eq.1.or.j.ne.1) k=k+1
0148              dphi(k)=dx(j)
0149              ddephi(k)=dh
0150              dsumph(k)=dsum+dz(j)
0151         enddo
0152         dsum=dsum+dz(ndim)
0153       enddo
0154       kmp=k
0155       return
0156       end
0157 
0158 C **********************************************************************
0159 C
0160 C     RADGAM_POL
0161 C
0162 C **********************************************************************
0163 
0164       subroutine radgam_pol(erad,yrad,xrad,genphi,
0165      +     ixytest,itagen,q2tr,anutr)
0166 
0167       include "double.inc"
0168       include "phiout.inc"
0169       include "radgen.inc"
0170       include "cmpcom.inc"
0171       include "sxycom.inc"
0172       include "tailcom.inc"
0173       include "amf2com.inc"
0174       include "density.inc"
0175 
0176       real*8 podinl
0177       real temp             ! Temp variable for a denom
0178 
0179 C     selection of mass dmj and mass index iadmj
0180       r1=rlu(0)
0181       if(ixytest.lt.0)then
0182         rr1=r1*(tbor+tine+tpro+tnuc)
0183 C      print *,tine,tpro,tnuc
0184         if (rr1.gt.(tine+tpro+tnuc)) then
0185           ita=0
0186         elseif(rr1.gt.(tpro+tnuc)) then
0187           ita=1
0188           sicurr=rr1-tpro-tnuc
0189           call conk2(dble(erad),dble(xrad),dble(yrad),1)
0190         elseif(rr1.gt.tnuc)then
0191           ita=3
0192           sicurr=rr1-tnuc
0193           call conk2(dble(erad),dble(xrad),dble(yrad),1)
0194         else
0195           ita=2
0196           sicurr=rr1
0197           call conk2(dble(erad),dble(xrad),dble(yrad),2)
0198         endif
0199       else
0200         xs=dble(xrad)
0201         ys=dble(yrad)
0202         ita=itagen
0203         if(ita.eq.1)sicurr=r1*tine
0204         if(ita.eq.2)sicurr=r1*tnuc
0205         if(ita.eq.3)sicurr=r1*tpro
0206         if(ita.eq.2)then
0207           call conk2(dble(erad),dble(xrad),dble(yrad),2)
0208         else
0209           call conk2(dble(erad),dble(xrad),dble(yrad),1)
0210 C           print *,an,xs,ys,ita,s,x
0211         endif
0212         isf1=1
0213         isf2=isf20
0214         isf3=1
0215 
0216 C      print *,an,xs,ys,erad,xrad,yrad
0217 C      sicurr=scgen
0218       endif
0219 
0220       if(ita.ne.0)then
0221 *
0222 *---     selection of theta angle
0223 *
0224 * maximum value of dsitkm can be smaller than 1 !
0225 * ---> use modified random number (as in old generator)
0226 * abr (03.01.02)
0227 *       print *,dsitkm(ndxtkm(ita),ita)
0228         r1=r1*dsitkm(ndxtkm(ita),ita)
0229         do ktk=1,ndxtkm(ita)
0230 C   print *,dsitkm(ktk,ita),r1,ndxtkm(ita)
0231           if (r1.le.dsitkm(ktk,ita)) then
0232             if (ktk.eq.1) then
0233               dtk1=r1/dsitkm(ktk,ita)
0234             else
0235               dtk1=(r1-dsitkm(ktk-1,ita))/
0236      +             (dsitkm(ktk,ita)-dsitkm(ktk-1,ita))
0237             endif
0238             dtk=dcvtkm(ktk,ita)+(dtk1-0.5d0)*ddetkm(ktk,ita)
0239             goto 30
0240           endif
0241         enddo
0242         dtk=dcvtkm(ndxtkm(ita),ita)
0243  30     continue
0244 * up to this point identical to old generator
0245 * exception: calculation of integrand
0246 *
0247         taout=(sx-sqly*cos(dtk))/ap2
0248 C       print *,taout,dtk,ita,ktk,r1
0249 C       print *,sx,sqly,ap2
0250         if(ita.eq.1)then
0251           taa=taout
0252           iphi=0
0253           call tails(taout,tm)
0254           rcal=ap*demin
0255           rmax=(w2-amc2)/(1.+taout)
0256           do krr=1,nrr
0257             ddeler(krr,ktk)=(rmax-rcal)/nrr
0258             drcurr(krr,ktk)=rcal+ddeler(krr,ktk)*(krr-0.5)
0259 C             print *,krr,ktk,drcurr(krr,ktk)
0260             if(krr.eq.1)then
0261              dsigmr(krr,ktk)=podinl(drcurr(krr,ktk))
0262             else
0263              dsigmr(krr,ktk)=dsigmr(krr-1,ktk)+podinl(drcurr(krr,ktk))
0264             endif
0265 C           print *,dsigmr(krr,ktk),drcurr(krr,ktk),taout
0266           enddo
0267           r2=rlu(0)
0268           dsircu=r2*dsigmr(nrr,ktk)
0269           do krr=1,nrr
0270             if (dsircu.le.dsigmr(krr,ktk)) then
0271               if(krr.eq.1)then
0272                 drr1=dsircu/dsigmr(krr,ktk)
0273               else
0274                 drr1=(dsircu-dsigmr(krr-1,ktk))
0275      +               /(dsigmr(krr,ktk)-dsigmr(krr-1,ktk))
0276               endif
0277               drr=drcurr(krr,ktk)+(drr1-0.5d0)*ddeler(krr,ktk)
0278               goto 20
0279             endif
0280           enddo
0281           drr=drcurr(nrr,ktk)
0282  20       continue
0283           dom=drr/ap
0284         else
0285           dom=(sx-y)/ap/(1.+taout)
0286         endif
0287         rrout=ap*dom
0288 
0289 * the following is again the same as in the old generator (abr 03.01.02)
0290 *
0291 *---selection of phi angle
0292 *
0293 *
0294 C  print *,rrout,taout
0295 C  pause
0296         iphi=1
0297         call phidst_pol
0298 *
0299         r3=rlu(0)
0300         dphpoi=r3*dsumph(kmp)
0301         do kph=2,kmp
0302           if (dphpoi.le.dsumph(kph)) then
0303 * avoid division by zero (abr)
0304 *           dphk1=(dphpoi-dsumph(kph-1))/
0305 *    +           (dsumph(kph)-dsumph(kph-1))
0306             temp = dsumph(kph) - dsumph(kph-1)
0307             if (temp .eq. 0) temp = 1e-20
0308             dphk1 = (dphpoi - dsumph(kph-1)) / temp
0309             dphk=dphi(kph)+(dphk1-1.0d0)*ddephi(kph)
0310             goto 40
0311           endif
0312         enddo
0313         dphk=dphi(kmp)
0314  40     continue
0315         r4=rlu(0)
0316         if (r4.gt.0.5) dphk=-dphk
0317         
0318 *
0319 *     radiative photon
0320 *
0321         deg=dom
0322         dthg=dtk
0323         dphig=dphk
0324         sigma_total=sngl(tbor+tine+tpro+tnuc)
0325         
0326         q2tr=y+rrout*taout
0327         anutr=sx/ap-dom
0328 C       write(6,*) 'ita,deg,dthg,dphig,q2tr,anutr,xtr',
0329 C     &              ita,deg,dthg,dphig,q2tr,anutr,
0330 C     &              q2tr/(2.*amp*anutr)
0331         
0332         dgpz=deg*dcos(dthg)
0333         dgpxy=deg*dsin(dthg)
0334         dgpx=dgpxy*dcos(dphig)
0335         dgpy=dgpxy*dsin(dphig)
0336 *     
0337 *---momentum components in the LAB-system:
0338 *---two rotations needed - first within the scattering plane
0339 *---around the y-axis and second around the new z-axis (beam
0340 *---direction) by Phi of the scattering plane
0341 *
0342         dgplx=-dgpz*dsts+dgpx*dcts
0343         dgply=dgpy
0344         dgplz=dgpz*dcts+dgpx*dsts
0345         dcphi=dcos(dble(genphi))
0346         dsphi=dsin(dble(genphi))
0347         dplabg(1)=dgplx*dcphi-dgply*dsphi
0348         dplabg(2)=dgplx*dsphi+dgply*dcphi
0349         dplabg(3)=dgplz
0350         
0351       else
0352 
0353         q2tr=2d0*amh*erad*xrad*yrad
0354         anutr=yrad*erad
0355 C         write(*,'(i5,5f8.3)')ita,xs,ys,q2tr,anutr,sigcor
0356 
0357         dplabg(1)=0.0d0
0358         dplabg(2)=0.0d0
0359         dplabg(3)=0.0d0
0360 
0361       endif
0362 
0363 * reset kinematics (abr)
0364 *     call conk2(dble(erad),dble(xrad),dble(yrad),1)
0365 * now done at beginning of sr radgen
0366 
0367       end
0368 
0369 C **********************************************************************
0370 C
0371 C     DFUFI_POL
0372 C
0373 C **********************************************************************
0374 
0375       double precision function dfufi_pol(dx)
0376 *     -----------------------------------
0377 *     needed for correction of hadronic distributions
0378       implicit real*8(a-h,o-z)
0379 
0380       include "cmpcom.inc"
0381       include "sxycom.inc"
0382       include "tailcom.inc"
0383       include "amf2com.inc"
0384       include "ppicom.inc"
0385       include "radgen.inc"
0386 
0387       phipoi=dx
0388       taa=taout
0389       call tails(taout,tm)
0390       dfufi_pol=an*alfa/pi*podinl(rrout)
0391       
0392       return
0393       end
0394       
0395 C **********************************************************************
0396 C
0397 C     RADQUAD_POL
0398 C
0399 C **********************************************************************
0400 
0401       subroutine radquad_pol(dfunct,dlow,dup,dx,dy,dz,dh,ndim)
0402 *     -------------------------------------------------
0403 *
0404       include "double.inc"
0405       dimension dx(31),dy(31),dz(31)
0406       external dfunct
0407 *
0408       dsum2=0.d0
0409       if (ndim.gt.1) then
0410         dh=(dup-dlow)/float(ndim-1)
0411         do i=1,ndim
0412           daux=dlow+dh*float(i-1)
0413           dx(i)=daux
0414           dy(i)=dfunct(daux)
0415         enddo
0416         do i=2,ndim
0417           dsum1=dsum2
0418           dsum2=dsum2+.5d0*(dx(i)-dx(i-1))*(dy(i)+dy(i-1))
0419           dz(i-1)=dsum1
0420         enddo
0421         dz(ndim)=dsum2
0422       elseif (ndim.eq.1) then
0423         dz(ndim)=dsum2
0424       endif
0425 *     
0426       return
0427       end
0428 
0429 C **********************************************************************
0430 C
0431 C     MPOLRAD
0432 C
0433 C **********************************************************************
0434 
0435       subroutine mpolrad(e1curr,yscurr,xscurr
0436      +     ,uncurr,plcurr,pncurr,itagen)
0437       implicit real*8(a-h,o-z)
0438 
0439       include "cmpcom.inc"
0440       include "sxycom.inc"
0441       include "ppicom.inc"
0442       include "tailcom.inc"
0443       include "radgen.inc"
0444       include "radgenkeys.inc"
0445       include "deltacom.inc"
0446       include "density.inc"
0447       include "pypars.inc"
0448 
0449       real e1curr,yscurr,xscurr,uncurr,plcurr,pncurr
0450 
0451       dimension tai(5),si(2,3),si0(2,3),tls(2,3,4)
0452 
0453       e1=e1curr
0454       xs=xscurr
0455       ys=yscurr
0456       pl=plcurr
0457       pn=-pncurr
0458       un=uncurr
0459 
0460       if (abs(pncurr).eq.2.) then
0461 c        negative polarization is twice as large as positive one
0462          qn=pncurr*1./2.
0463          if(qn.lt.0) qn=2.*qn
0464       else
0465          qn=0.
0466       endif
0467 
0468       call conk2(e1,xs,ys,1)
0469 c
0470 c
0471 c delta is factorizing part of virtual and real leptonic bremsstrahlung
0472 c
0473       call deltas(delta,delinf,tr)
0474 
0475       il=1
0476       in=1
0477 
0478       si(il,in)=0d0
0479       si0(il,in)=0d0
0480       tls(il,in,1)=0d0
0481       tls(il,in,2)=0d0
0482       tls(il,in,3)=0d0
0483       tls(il,in,4)=0d0
0484       
0485       isf1=1
0486       isf2=isf20
0487       isf3=1
0488 
0489       do ii=1,4
0490         tai(ii)=0d0
0491       enddo
0492       sib=0d0
0493       sia=0d0
0494 
0495       if(itagen.eq.-1)then
0496        ita1=1
0497        ita4=4
0498       else
0499        ita1=itagen
0500        ita4=itagen
0501       endif
0502 
0503       do 30 ita=ita1,ita4
0504 c
0505 c     sib is born cross section with polarized initial
0506 c     lepton and proton
0507 c     sia is contribution of anomalous magnetic moment.
0508 c
0509 * for photoproduction set elastic and quasielastic tails to zero
0510 * abr (19.02.04)
0511       if (ntx.eq.ntpho .and. (ita.eq.2.or.ita.eq.3)) then
0512         tai(2)=0.0d0
0513         tai(3)=0.0d0
0514         goto 30
0515       endif
0516 
0517 C ECA to use Antjes trick to kill the elastic and quasielastic tails
0518 C always if we run pythia with radiative corrections
0519 c      if (MSTP(199).eq.1) then
0520 c        tai(2)=0.0d0
0521 c        tai(3)=0.0d0
0522 c        goto 30
0523 c      endif
0524 
0525         if(ita.eq.2.and.kill_elas_res.eq.1)goto30
0526         if(ita.eq.3.and.ire.le.1)then
0527           tai(3)=0d0
0528 C           write(9,'('' tai = 0.0  '')')
0529           goto 30
0530         end if
0531         if(ita.eq.1)then
0532           call bornin(sib,sia)
0533         endif
0534 c
0535 c     tai(1),tai(2),tai(3) are contributions of radiative tails:
0536 c     1 - inelastic
0537 c     2 - elastic
0538 c     3 - quasielastic
0539 c
0540         if(ita.eq.2) call conk2(e1,xs,ys,ita)
0541 
0542         if(kill_elas_res.ne.2) call qqt(sib,tai(ita))
0543 
0544 * what is this ??? (abr 17.10.01)
0545 * shouldnot this if statement be not equal here ???
0546 * abr   if(ita.eq.2) call conk2(e1,xs,ys,1)
0547 *       if(ita.ne.2) call conk2(e1,xs,ys,1)
0548 * no, purpose is to reset kinematics, I think (abr 03.01.02)
0549 * thus the original code is correct
0550         if(ita.eq.2) call conk2(e1,xs,ys,1)
0551 
0552  30   continue
0553 C       extai1=exp(alfa/pi*delinf)
0554 C       extai2=((sx-y/rtara)**2/s/(s-y/rtara))**tr
0555 C       extai3=((sx-y)**2/s/(s-y))**tr
0556 
0557       extai1=1.
0558       extai2=1.
0559       extai3=1.
0560       delinf=0.
0561 
0562 C       sig=sib*extai1*(1.+alfa/pi*(delta-delinf))+sia
0563 C      +       +tai(4)+tai(5)
0564 cilyichev only one-loop without multy-soft photon emission
0565 c      sig=sib*redfac*(1.+vertex+vac+small)+sia
0566       sig=sib*(1.+log(redfac)+vertex+vac+small)+sia
0567      +     +tai(4)
0568      +     +tai(1)+(tai(2)*extai2+tai(3)*extai3)/rtara
0569 
0570       si(il,in)=si(il,in)+sig
0571       si0(il,in)=si0(il,in)+sib
0572       tls(il,in,1)=tls(il,in,1)+tai(1)
0573       tls(il,in,2)=tls(il,in,2)+tai(2)*extai2/rtara
0574       tls(il,in,3)=tls(il,in,3)+tai(3)*extai3/rtara
0575 C       tls(il,in,4)=sib*extai1*(1.+alfa/pi*(delta-delinf))+sia
0576 C      +     +tai(4)+tai(5)
0577       tls(il,in,4)=sib*redfac*(1.+vertex+vac+small)+sia
0578      +     +tai(4)
0579 
0580 C      write(*,'(1x,8g11.4)')xs,ys,si0(il,in),si(il,in)
0581 C      +  ,tls(il,in,1),tls(il,in,2),tls(il,in,3)
0582 C      +  ,tls(il,in,4)
0583       if(itagen.eq.-1.or.itagen.eq.1)then
0584 C        if(abs(tine-tls(il,in,1))/tine.gt.5d-2)
0585 C      +     write(*,*)' tine',tine,tls(il,in,1),itagen
0586         tine=tls(il,in,1)
0587       endif
0588       if(itagen.eq.-1.or.itagen.eq.2)then
0589 C        if(abs(tnuc-tls(il,in,2))/tnuc.gt.5d-2)
0590 C      +     write(*,*)' tnuc',tnuc,tls(il,in,2),itagen
0591         tnuc=tls(il,in,2)
0592       endif
0593       if(itagen.eq.-1.or.itagen.eq.3)then
0594 C      if(abs(tpro-tls(il,in,3))/tpro.gt.5d-2)
0595 C      +     write(*,*)' tpro',tnuc,tls(il,in,3),itagen
0596         tpro=tls(il,in,3)
0597       endif
0598       if(itagen.eq.-1)then
0599         tbor=tls(il,in,4)
0600         sigrad=si(il,in)
0601         sig1g=si0(il,in)
0602       endif
0603 C #ifdef TEST
0604 C       write(*,'(1x,6g11.4)')sig1g,sigrad,tine,tnuc,tpro,tbor
0605 C #endif
0606       sigcor=sigrad/sig1g
0607       if(kill_elas_res.eq.2)sigcor=1.
0608       
0609 C #ifdef TEST
0610 C       print *,' kill_elas_res=',kill_elas_res,amc2
0611 C       print *,' sig1g=',sig1g,sib
0612 C       print *,' sigrad=',sigrad,sig
0613 C       print *,' tine =',tine
0614 C       print *,' tai(4) =',tai(4)
0615 C       print *,' tnuc =',tnuc
0616 C       print *,' tpro =',tpro
0617 C       print *,' tbor =',tbor
0618 C       print *,' sig1g=',sig1g
0619 C       print *,' sigcor=',sigcor
0620 C       print *,' delta=',alfa/pi*delta
0621 C       print *,' vac  =',vac
0622 C       print *,' vertex=',vertex
0623 C       print *,' small=',small
0624 C       print *,' tai(4)= ',sib*log(redfac)
0625 C       sigt1=sib*alfa/pi*delta+tai(5)
0626 C       sigt2=sib*(log(redfac)+vertex+vac+small)
0627 C       sigt3=sib*alfa/pi*(delta+delta5)
0628 C       print *,' sigt=',sigt1,sigt2,sigt3
0629 C #endif
0630 *     print *,' kill_elas_res=',kill_elas_res,amc2
0631 *     print *,' sig1g=',sig1g,sib
0632 *     print *,' sigrad=',sigrad,sig
0633 *     print *,' tine =',tine
0634 *     print *,' tai(4) =',tai(4)
0635 *     print *,' tnuc =',tnuc
0636 *     print *,' tpro =',tpro
0637 *     print *,' tbor =',tbor
0638 *     print *,' sig1g=',sig1g
0639 *     print *,' sigcor=',sigcor
0640 *     print *,' delta=',alfa/pi*delta
0641 *     print *,' vac  =',vac
0642 *     print *,' vertex=',vertex
0643 *     print *,' small=',small
0644 *     print *,' redfac= ',redfac
0645 
0646       end
0647 
0648 C **********************************************************************
0649 C
0650 C     CONK2
0651 C
0652 C **********************************************************************
0653 
0654       subroutine conk2(e1,xs,ys,iittaa)
0655       implicit real*8(a-h,o-z)
0656 
0657       include "cmpcom.inc"
0658 
0659       if(iittaa.eq.2)then
0660         amp=amt
0661       else
0662         amp=amh
0663       endif
0664 C      print *,amp,amh,iittaa
0665       call conkin(e1,xs*amh/amp,ys)
0666 
0667       return
0668       end
0669 
0670 C **********************************************************************
0671 C
0672 C     CONKIN
0673 C
0674 C **********************************************************************
0675 
0676       subroutine conkin(e1,xss,yss)
0677 
0678 c set of kinematical constants
0679       implicit real*8(a-h,o-z)
0680 
0681       include "cmpcom.inc"
0682       include "polcom.inc"
0683       include "sxycom.inc"
0684       include "ppicom.inc"
0685       include "radgen.inc"
0686 cilyichev additional include
0687       include "mc_set.inc"
0688 
0689       ap=2.*amp
0690       amp2=amp**2
0691       ap2=2.*amp**2
0692       s=ap*e1
0693       x=s*(1-yss)
0694       sx=s-x
0695       sxp=s+x
0696       y=s*xss*yss
0697       ym=y+al2
0698       tpl=s**2+x**2
0699       tmi=s**2-x**2
0700       w2=amp2+s-y-x
0701       als=s*s-al2*ap2
0702       alx=x*x-al2*ap2
0703       alm=y*y+4.*aml2*y
0704       aly=sx**2+4.*amp2*y
0705       sqls=dsqrt(als)
0706       sqlx=dsqrt(alx)
0707       sqly=dsqrt(aly)
0708       sqlm=dsqrt(alm)
0709       allm=dlog((sqlm+y)/(sqlm-y))/sqlm
0710 
0711       coe=xs/e1/1d3
0712 
0713       axy=pi*(s-x) * coe
0714       an=2.*alfa**2/sqls*axy*barn*amh/amp
0715       tamin=(sx-sqly)/ap2
0716       tamax=(sx+sqly)/ap2
0717 
0718       dcts=(s*sx + ap2*y)/sqly/sqls
0719 * avoid values above 1.0 (for x values below 1e-07, abr 18.2.04)
0720 *     dsts=sin( acos(dcts) )
0721       dsts=sin( acos(MIN(MAX(dcts,-1.0),+1.0)) )
0722 
0723       as=s/2./aml/sqls
0724       bs=0.
0725       cs=-aml/sqls
0726 cilyichev instead of
0727 c        ae=amp/sqls
0728 c        be=0.
0729 c        ce=-s/ap/sqls
0730 c we need
0731       if ((mcSet_PTarget(1:1).eq.'L').or.
0732      +    (mcSet_PTarget(1:2).eq.'DT')) then
0733          ae=amp/sqls
0734          be=0.
0735          ce=-s/ap/sqls
0736       elseif(mcSet_PTarget(1:1).eq.'T')then
0737          sqn=dsqrt(s*x*y-aly*aml2-amp2*y*y)
0738          ae=(-s*x+ap2*ym)/sqls/sqn/2.
0739          be=sqls/sqn/2.
0740          ce=-(s*y+al2*sx)/sqls/sqn/2.
0741       endif
0742 
0743       apq=-y*(ae-be)+ce*sx
0744       apn=(y+4.*aml2)*(ae+be)+ce*sxp
0745       dk2ks=as*ym+al2*bs+cs*x
0746       dksp1=as*s+bs*x+cs*ap2
0747       dapks=2.*(al2*(as*ae+bs*be)+ap2*cs*ce+ym*(as*be+bs*ae)
0748      +     +s*(as*ce+cs*ae)+x*(bs*ce+cs*be))
0749 
0750       return
0751       end
0752 
0753 C **********************************************************************
0754 C
0755 C     BORNIN
0756 C
0757 C **********************************************************************
0758 
0759       subroutine bornin(sibor,siamm)
0760 c
0761 c     sibor is born cross section with polarized initial
0762 c     lepton and polarized target
0763 c     siamm is contribution of anomalous magnetic moment.
0764 C     the cross section is calculated as dsigma/dlogQ2dnu to fit
0765 C     to the HERMES-disng generation
0766 
0767 c
0768       implicit real*8(a-h,o-z)
0769 
0770       include "cmpcom.inc"
0771       include "polcom.inc"
0772       include "sxycom.inc"
0773       include "ppicom.inc"
0774       include "tailcom.inc"
0775 
0776       dimension sfm0(8),tm(8)
0777 
0778       call strf(0d0,0d0,sfm0)
0779       tm(1)=-(2.*aml2-y)
0780       tm(2)=(-(amp2*y-s*x))/(2.*amp2)
0781       tm(3)=(2.*(apq*dk2ks-dapks*y)*aml)/amp
0782       tm(4)=apq/amp*(-(dk2ks*sx-2.*dksp1*y)*aml)/amp2
0783       tm(7)=(-(4.*aml2+3.*apn**2-3.*apq**2+y))/2.
0784       tm(8)=apq/amp*(-3.*(apn*sxp-apq*sx))/(2.*ap)
0785       ek=(3.*apq**2-y)/amp2
0786       tm(5)=-ek*tm(1)
0787       tm(6)=-ek*tm(2)
0788       ssum=0.
0789       do 1 isf=isf1,isf2,isf3
0790          ppol=1.
0791          if(isf.eq.3.or.isf.eq.4)ppol=pl*pn
0792 C        if(isf.eq.3.or.isf.eq.4)ppol=-pn
0793          if(isf.ge.5)ppol=qn/6
0794         ssum=ssum+tm(isf)*sfm0(isf)*ppol
0795     1 continue
0796       sibor=ssum*an/y/y*2.
0797 c     
0798 c     formula (4) of kukhto and shumeiko paper
0799 c
0800 cc    res1=amp*ww1*(y+4.*aml2)-ww2*(s+x)**2/4./amp
0801 cc    siamm=alfa/pi*al2*allm*(sibor+an*res1/y**2)
0802       siamm=0.
0803 
0804       return
0805       end
0806 
0807 C **********************************************************************
0808 C
0809 C     DELTAS
0810 C
0811 C **********************************************************************
0812 
0813       subroutine deltas(delta,delinf,tr)
0814 c
0815 c delta is factorizing part of virtual and real leptonic bremsstrahlung
0816 c
0817       implicit real*8(a-h,o-z)
0818 
0819       include "cmpcom.inc"
0820       include "sxycom.inc"
0821       include "ppicom.inc"
0822       include "radgen.inc"
0823       include "deltacom.inc"
0824 
0825 c
0826 c    am2 : squared masses of charge leptons
0827 c
0828       dimension am2(3)
0829       data am2/.26110d-6,.111637d-1,3.18301d0/
0830 
0831       del1=-ym*(alm*allm**2/2.+2.*fspen(2d0*sqlm/(y+sqlm))-pi2/2.)/sqlm
0832       del2=(3.*y/2.+4.*aml2)*allm-2.
0833 
0834       suml=0.
0835       do 10 i=1,3
0836         a2=2.*am2(i)
0837         sqlmi=dsqrt(y*y+2.*a2*y)
0838         allmi=dlog((sqlmi+y)/(sqlmi-y))/sqlmi
0839         suml=suml+2.*(y+a2)*allmi/3.-10./9.+4.*a2*(1.-a2*allmi)/3./y
0840  10   continue
0841       if(y.lt.1.d0)then
0842         aaa = -1.345d-9
0843         bbb = -2.302d-3
0844         ccc = 4.091
0845       elseif(y.lt.64d0)then
0846         aaa = -1.512d-3
0847         bbb =  -2.822d-3
0848         ccc = 1.218
0849       else
0850         aaa = -1.1344d-3
0851         bbb = -3.0680d-3
0852         ccc = 9.9992d-1
0853       endif
0854       sumh = -(aaa+bbb*log(1.+ccc*y)) *2*pi/alfa
0855       sum=suml+sumh
0856 C       print *,' vacl_my_old=',alfa/pi*suml
0857 C       print *,' vach_my_old=',alfa/pi*sumh
0858 C       print *,' vac_my_old=',alfa/pi*sum
0859       sum=vacpol(y)
0860 C       print *,' vac_my_new=',alfa/pi*sum
0861 
0862       aj0=2.*(ym*allm-1.)
0863 c...elke okay the dlog((w2-amc2) can become negative if w2 is smaller than
0864 C   amc2/1.151857d0/
0865 c     write(*,*)w2,aml,amc2,aj0
0866       deltai=aj0*dlog((w2-amc2)/aml/dsqrt(w2))
0867 c     write(*,*)"deltai=",deltai
0868 
0869       ss=x+y
0870       xx=s-y
0871       alss=ss**2-2.*w2*al2
0872       alxx=xx**2-2.*w2*al2
0873       sqlss=dsqrt(alss)
0874       sqlxx=dsqrt(alxx)
0875       allss=dlog((sqlss+ss)/(-sqlss+ss))/sqlss
0876       allxx=dlog((sqlxx+xx)/(-sqlxx+xx))/sqlxx
0877       dlm=dlog(y/aml2)
0878       sfpr=dlm**2/2.-dlm*dlog(ss*xx/aml2/w2)
0879      +     -(dlog(ss/xx))**2/2.+fspen((s*x-y*amp2)/ss/xx)-pi2/3.
0880       delta0=(ss*allss+xx*allxx)/2.+sfpr
0881       delta=deltai+delta0+del1+del2+sum
0882       delinf=(dlm-1.)*dlog((w2-amc2)**2/ss/xx)
0883       tr=alfa/pi*(dlm-1.)
0884 
0885       vac=alfa/pi*sum
0886       vertex=alfa/pi*del2
0887       small_old=alfa/pi*(pi2/6.-fspen(1.-amp2*y/s/x)
0888      +     + fspen(1.-s/x) + fspen(1.-x/s))
0889       small=alfa/pi*(-pi2/6.-log(ss*xx/s/x)*log(sx/y)
0890      +     +log(ss/s)*log(xx/x)
0891      +     -1.5*log(s/x)**2+2.*log(xx/x)*log(s/y)+2.*log(ss/s)*log(x/y)
0892      +     -2.*fspen(-y/x)-2.*fspen(y/s)+2.*fspen(y/sx)+fspen(s*x/ss/xx)
0893      +     -fspen(s*y/ss/sx)-fspen(x*y/sx/xx))
0894 C       print *,' small_old=',small_old
0895 C       print *,' small_new=',small
0896         redfac=exp(-alfa/pi*(dlm-1.)*log(s*x/(4.*amp2*demin**2) ))
0897         delta5=-(dlm-1.)*log((w2-amc2)**2*s*x/(4.*amp2*demin**2*ss*xx))
0898      +       -log(xx/x)*log(amp2*y/s**2)-log(ss/s)*log(amp2*y/x**2)
0899      +       -2.*fspen(-y/x)-2.*fspen(y/s)+2.*fspen(-tamin)
0900      +       +2.*fspen(-tamax)
0901      +       -fspen((-y-s*tamin)/xx)
0902      +       -fspen((-y-s*tamax)/xx)
0903      +       -fspen(( y-x*tamin)/ss)
0904      +       -fspen(( y-x*tamax)/ss)
0905 
0906       return
0907       end
0908 
0909 C **********************************************************************
0910 C
0911 C     VACPOL
0912 C
0913 C **********************************************************************
0914 
0915       double precision function vacpol(y)
0916 
0917       implicit real*8(a-h,o-z)
0918 
0919       include "ppicom.inc"
0920 
0921 c
0922 c    am2 : squared masses of charge leptons
0923 c
0924       dimension am2(3),a(5),b(5),c(5)
0925       data am2/.26110d-6,.111637d-1,3.18301d0/
0926       data a/0d0,0d0,0d0,1.2227d-3,1.64178d-3/
0927       data b/2.2877d-3,2.51507d-3,2.79328d-3,2.96694d-3,2.92051d-3/
0928       data c/4.08041425d0,3.09624477d0,2.07463133d0,1d0,1d0/
0929 
0930       suml=0.
0931       do 10 i=1,3
0932         a2=2.*am2(i)
0933         sqlmi=dsqrt(y*y+2.*a2*y)
0934         allmi=dlog((sqlmi+y)/(sqlmi-y))/sqlmi
0935         suml=suml+2.*(y+a2)*allmi/3.-10./9.+4.*a2*(1.-a2*allmi)/3./y
0936  10   continue
0937 
0938       if(y .lt. 4d0)then
0939        k=1
0940       elseif(y .lt. 16d0)then
0941        k=2
0942       elseif(y .lt. 100d0)then
0943        k=3
0944       elseif(y .lt. 8317.44d0)then
0945        k=4
0946       elseif(y .ge. 8317.44d0)then
0947        k=5
0948       else
0949        stop ' Y<0 in VACPOL'
0950       endif
0951 
0952       sumh = (a(k)+b(k)*log(1.+c(k)*y)) *2.*pi/alfa
0953       vacpol=suml+sumh
0954 
0955       end
0956 
0957 C **********************************************************************
0958 C
0959 C     FSPENS
0960 C
0961 C **********************************************************************
0962 
0963       double precision function fspens(x)
0964 c
0965 c    spence function
0966 c
0967       implicit real*8(a-h,o-z)
0968 
0969       f=0.d0
0970       a=1.d0
0971       an=0.d0
0972       tch=1.d-16
0973  1    an=an+1.d0
0974       a=a*x
0975       b=a/an**2
0976       f=f+b
0977       if(b-tch)2,2,1
0978  2    fspens=f
0979 
0980       return
0981       end
0982 
0983 C **********************************************************************
0984 C
0985 C     FSPEN
0986 C
0987 C **********************************************************************
0988 
0989       double precision function fspen(x)
0990 
0991       implicit real*8(a-h,o-z)
0992 
0993       data f1/1.644934d0/
0994 
0995       if(x)8,1,1
0996  1    if(x-.5d0)2,2,3
0997  2    fspen=fspens(x)
0998       return
0999  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         print *, nny, nty, nt, ys, y(nny)
2180         if (ys.gt.y(nny)) goto 10
2181       enddo
2182  5    print *,' ys=',ys,' out of the region'
2183       stop
2184  10   if (nny.eq.nt) goto 5
2185 
2186       if(ixytest.eq.1)then
2187         stop ' ixytest.eq.1 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 **********************************************************************
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!