Back to home page

EIC code displayed by LXR

 
 

    


Warning, file /LQGENEP/lqgenep.f was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

0001 *-----------------
0002 * File: lqgenep.f
0003 *-----------------
0004 *
0005            subroutine LQGENEP(Nevt,flag)
0006 C------------------------------------------
0007 C...Main program for leptoquark generation 
0008 C...in electron-proton scattering
0009 C------------------------------------------
0010 C...All real arithmetic in double precision.
0011       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0012       Integer flag, itau,id,ii,ij
0013 
0014 C...LQGENEP run setup parameters
0015       double precision BEAMPAR,LQGPAR3,
0016      >  ptnt,phnt,ptt,pht,pth,phh,
0017      >  ptx,pty,ptz,phx,phy,phz,
0018      > ptid, phid,ppid,pxid,pyid,pzid
0019 
0020       CHARACTER*256 out_file
0021       integer LQGPAR1,LQGPAR2
0022       COMMON/LQGDATA/BEAMPAR(3),LQGPAR1(10),LQGPAR2(10),LQGPAR3(20),
0023      > out_file
0024 
0025 C...LQGENEP event informations
0026        double precision LQGKPAR,LQGDDX
0027        integer LQGPID
0028        COMMON/LQGEVT/LQGKPAR(3),LQGDDX(3),LQGPID(3)
0029 
0030 C...Pythia declarations. 
0031 C...Three Pythia functions return integers, so need declaring.
0032       INTEGER PYK,PYCHGE,PYCOMP
0033 C...Parameter statement to help give large particle numbers
0034 C...(left- and righthanded SUSY, excited fermions).
0035       PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KEXCIT=4000000)
0036 
0037 C...EXTERNAL statement links PYDATA on most machines.
0038       EXTERNAL PYDATA
0039 
0040 *
0041 C...Pythia Commonblocks.
0042 C...The event record.
0043       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0044 C...Parameters.
0045       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0046 C...Particle properties + some flavour parameters.
0047       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0048 C...Decay information.
0049       COMMON/PYDAT3/MDCY(500,3),MDME(4000,2),BRAT(4000),KFDP(4000,5)
0050 C...Selection of hard scattering subprocesses.
0051       COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
0052 C...Parameters. 
0053       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0054 C...Process information.
0055       COMMON/PYINT1/MINT(400),VINT(400)
0056       COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
0057 C...Supersymmetry parameters.
0058       COMMON/PYMSSM/IMSS(0:99),RMSS(0:99)
0059       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYDAT3/,/PYSUBS/,/PYPARS/,
0060      >/PYINT2/,/PYMSSM/
0061 C...Local Array.
0062       DIMENSION NCHN(12),QVEC(4)
0063       DATA NCHN/12*0/
0064 
0065 C...Internal used common
0066 C# LQGpar1.inc #
0067       integer echar
0068       double precision ebeam,pbeam
0069       common /LQGbeam/ebeam,pbeam,echar
0070 
0071 C# LQGpdfC.inc #
0072       character*20 parm(20)
0073       double precision pdfsf(20)
0074       common /LQGpdfC/ pdfsf 
0075       
0076 C# LQGKinC.inc #
0077       double precision xmax,xmin,ymax,ymin,zmax,zmin,Q2min
0078       common /LQGKinC/ xmax,xmin,ymax,ymin,zmax,zmin,Q2min
0079 
0080 C# LQGproc.inc #
0081       double precision Mlq,G1,G2
0082       Integer LQtype,l_o,q_i,q_j
0083       common /LQGproc/ Mlq,G1,G2,LQtype,l_o,q_i,q_j 
0084 
0085 C# LQGKinV.inc # 
0086       double precision S,Srad,x,y,z,Q2
0087       common /LQGKinV/ S,Srad,x,y,z,Q2
0088 
0089 C# LQGout.inc #
0090       double precision DXSEC(3),pvalence
0091       integer q_o,q_s,genproc,genprtyp,sch,uch,gproc(8)
0092       common /LQGout/ DXSEC,pvalence,q_o,q_s,genproc,genprtyp,
0093      >sch,uch,gproc
0094 
0095 
0096 C...processes
0097       integer ibea
0098       Character*9 chbea(2)
0099       Character*12 chprod(2,3)
0100       data chbea /'e+ qi -> ','e- qi -> '/
0101       data chprod /' -> e+ qj',' -> e- qj',
0102      >             ' -> mu+ qj',' -> mu- qj',
0103      >             ' -> tau+ qj',' -> tau- qj'/
0104 C...LQ type
0105       Character*7 LQCHA(14)
0106       DATA LQCHA /'S_0L','S_0R','~S_0R','S_1L',
0107      >            'V_1/2L','V_1/2R','~V_1/2L',
0108      >            'V_0L','V_0R','~V_0R','V_1L',
0109      >            'S_1/2L','S_1/2R','~S_1/2L'/                        
0110 
0111 C...Local declarations
0112       Real pxtot,pytot,pztot,etot,ecmtot,
0113      >     pxsum,pysum,pzsum,esum,ecmsum,
0114      >     pxlow,pxhi,pylow,pyhi,pzlow,pzhi,elow,ehi,ecmlow,ecmhi 
0115 
0116 C...ASCII Output File
0117       INTEGER asciiN
0118       PARAMETER (asciiN=29)
0119 !      CHARACTER*256 out_file
0120 
0121 C-----------------------------------------------------------------
0122 
0123       Integer Nwds_HBOOK
0124       Parameter (Nwds_HBOOK=100000)
0125       Real          HMEM
0126       Common /PAWC/ HMEM(Nwds_HBOOK)
0127 
0128 C Give output file name
0129 
0130 !      out_file='./TestOut.txt'
0131 C      out_file='./outdir/TestOut.txt'
0132 
0133 C...FLAG=0 -> First section: inizialization  
0134       If(flag.eq.0)then
0135 
0136 C.. LQGENEP banner
0137        call LQGBAN
0138 
0139 C.. Hbook inizialization       
0140        if(LQGPAR1(3).gt.0)Call HLIMIT(Nwds_HBOOK)
0141        if(LQGPAR1(3).lt.0)Call HLIMIT(-Nwds_HBOOK)
0142 
0143 C...beams properties.
0144        echar=beampar(1)
0145        Ebeam=beampar(2)
0146        Pbeam=beampar(3)
0147        S=4.d0*Ebeam*Pbeam
0148 
0149 C...LQ properties
0150        MLQ=LQGPAR3(1)
0151        G1=LQGPAR3(2)
0152        G2=LQGPAR3(3)
0153        LQTYPE=LQGPAR2(1)       
0154 
0155 C... outcoming lepton
0156        l_o=LQGPAR2(4)
0157 
0158 C... incoming and outcoming quark generation
0159        q_i=LQGPAR2(2)       
0160        q_j=LQGPAR2(3)
0161 
0162 C... kinematic ranges
0163        xmin=LQGPAR3(4)              
0164        xmax=LQGPAR3(5)              
0165        ymin=LQGPAR3(6)              
0166        ymax=LQGPAR3(7) 
0167        Zmin=0.d0
0168        Zmax=1.d0
0169        Q2min=LQGPAR3(8)             
0170 
0171 C... print LQGENEP generation run settings
0172        call LQGPRI1
0173 
0174 C... structure function
0175        parm(1)='NPTYPE'
0176        parm(2)='NGROUP'
0177        parm(3)='NSET'
0178        pdfsf(1)= LQGPAR3(9)
0179        pdfsf(2)= LQGPAR3(10) 
0180        pdfsf(3)= LQGPAR3(11)
0181        call PDFSET(PARM,pdfsf)
0182 
0183 C... Pythia initialization 
0184        P(1,1)=0D0
0185        P(1,2)=0D0
0186        P(1,3)=-Ebeam
0187        P(2,1)=0D0
0188        P(2,2)=0D0
0189        P(2,3)=Pbeam
0190 
0191 C...Evaluate limits for total momentum and energy histos 
0192        if(LQGPAR1(3).gt.0)then
0193         pxtot=sngl(P(1,1)+P(2,1))
0194         pytot=sngl(P(1,2)+P(2,2))
0195         pztot=sngl(P(1,3)+P(2,3))
0196         etot=sngl(dabs(P(1,3))+dabs(P(2,3)))
0197         ecmtot=sqrt(etot*etot-(pxtot*pxtot+pytot*pytot+pztot*pztot))
0198         if(pxtot.gt.0)then
0199          pxlow=pxtot-0.01*pxtot
0200          pxhi=pxtot+0.01*pxtot
0201         else
0202          pxlow=pxtot-1.
0203          pxhi=pxtot+1.
0204         endif
0205         if(pytot.gt.0)then
0206          pylow=pytot-0.01*pytot
0207          pyhi=pytot+0.01*pytot
0208         else
0209          pylow=pytot-1.
0210          pyhi=pytot+1.
0211         endif
0212         if(pztot.gt.0)then
0213          pzlow=pztot-0.01*pztot
0214          pzhi=pztot+0.01*pztot
0215         else
0216          pzlow=pztot-1.
0217          pzhi=pztot+1.
0218         endif
0219         if(etot.gt.0)then
0220          elow=etot-0.01*etot
0221          ehi=etot+0.01*etot
0222         else
0223          elow=etot-1.
0224          ehi=etot+1.
0225         endif                        
0226         if(ecmtot.gt.0)then
0227          ecmlow=ecmtot-0.01*ecmtot
0228          ecmhi=ecmtot+0.01*ecmtot
0229         else
0230          ecmlow=ecmtot-1.
0231          ecmhi=ecmtot+1.
0232         endif        
0233        endif
0234 C...Initialize Pythia
0235        if(LQGPAR1(5).eq.0)then
0236         Isub=401
0237        else
0238         Isub=LQGPAR1(5)
0239        endif
0240        sch=0.d0
0241        uch=0.d0
0242        call vzero(gproc,8)
0243        sigmax=LQGPAR3(12)
0244        if(echar.gt.0)then
0245         ibea=1
0246        else
0247         ibea=2
0248        endif
0249 
0250        
0251        CALL PYUPIN(ISUB,
0252      >      CHBEA(ibea)//LQCHA(LQTYPE)//CHPROD(ibea,l_o),sigmax)
0253        MSEL=0
0254        MSUB(ISUB)=1
0255 *
0256        if(beampar(1).GT.0)then
0257         CALL PYINIT('USER','e+','p',0D0)
0258        else
0259         CALL PYINIT('USER','e-','p',0D0)
0260        endif
0261 
0262 
0263        if(LQGPAR1(3).ne.0)then
0264 C...Book histos.
0265         call hropen(69,'lqgenep','lqgenep.histo','N',1024,ierr) 
0266         CALL hbook1(1000,'x gen',50,sngl(xmin),sngl(xmax),0.)
0267         CALL hbook1(1001,'x gen s-ch.',50,sngl(xmin),sngl(xmax),0.)
0268         CALL hbook1(1002,'x gen u-ch',50,sngl(xmin),sngl(xmax),0.)
0269         CALL hbook1(2000,'y gen',50,sngl(ymin),sngl(ymax),0.) 
0270         CALL hbook1(2001,'y gen s-ch',50,sngl(ymin),sngl(ymax),0.) 
0271         CALL hbook1(2002,'y gen u-ch',50,sngl(ymin),sngl(ymax),0.) 
0272         CALL hbook1(3000,'Q2 gen',50,0.,6.,0.) 
0273         call hbook1(5001,'sum px',100,pxlow,pxhi,0.)
0274         call hbook1(5002,'sum py',100,pylow,pyhi,0.) 
0275         call hbook1(5003,'sum pz',100,pzlow,pzhi,0.)
0276        call hbook1(5004,'sum e',100,elow,ehi,0.) 
0277         call hbook1(5000,'center of mass energy',
0278      >  100,ecmlow,ecmhi,0.) 
0279        endif
0280        write(6,*)
0281       endif
0282 
0283 C...Open File
0284       OPEN(asciiN, file=out_file)
0285 
0286 C      WRITE(29,*)'============================================'
0287       Integer HeaderTest = 0
0288       if(HeaderTest .LE. 1) then
0289       WRITE(29,*)' Pythia Output File'
0290       WRITE(29,*)'============================================'
0291       WRITE(29,*)'I, ievent, genevent, subprocess, nucleon,
0292      & targetparton, xtargparton, beamparton, xbeamparton,
0293      & thetabeamprtn, truey, trueQ2, truex, trueW2, trueNu,
0294      & leptonphi, s_hat, t_hat, u_hat, pt2_hat, Q2_hat, F2, F1, R,
0295      &  sigma_rad, SigRadCor, EBrems, photonflux, nrTracks'
0296       WRITE(29,*)'============================================'
0297 C      WRITE(29,30) 'I', 'KS', 'KF(ID)', 'ORIG'
0298 C30    FORMAT(A5, A5, A7, A5)
0299 C      WRITE(29,30) 'I', 'KS', 'KF(ID)', 'ORIG', 'px'
0300 C     &           , 'py', 'pz', 'E', 'm'
0301 C30    FORMAT(A5, A5, A7, A5, A10, A10, A10, A10, A10)
0302       WRITE(29,30) 'I', 'KS', 'KF(ID)', 'ORIG', 'D1', 'D2', 'px'
0303      &           , 'py', 'pz', 'E', 'm', 'vx', 'vy', 'vz'
0304 30    FORMAT(A5, A5, A7, A5, A5, A5, A10, A10, A10, A10, A10, 
0305      &           A10, A10, A10)
0306       WRITE(29,*)'============================================'
0307       endif
0308       HeaderTest = 999
0309 C-----------------------------------------------------------------
0310 C...FLAG=1 -> Second section: event generation  
0311 
0312       if(flag.eq.1)Then
0313        CALL PYEVNT
0314 c       print*,"The no. of event is",LQGPAR1(4)
0315        CALL PYHEPC(1)
0316 
0317 C...s-u channel
0318        if(genproc.eq.1)sch=sch+1     
0319        if(genproc.eq.2)uch=uch+1     
0320 
0321 C...process type
0322       gproc(genprtyp)=gproc(genprtyp)+1
0323 
0324 C...Fill event informations common
0325        LQGKPAR(1)=X
0326        LQGKPAR(2)=Y
0327        LQGKPAR(3)=Q2
0328        LQGDDX(1)=(DXSEC(2)+DXSEC(3))*1.d-9
0329        LQGDDX(2)=DXSEC(3)*1.d-9
0330        LQGDDX(3)=DXSEC(2)*1.d-9
0331        LQGPID(1)=q_s
0332        LQGPID(2)=q_o
0333        if(genproc.eq.1)then
0334         LQGPID(3)=1
0335        elseif(genproc.eq.2)then
0336         LQGPID(3)=2
0337        endif
0338 
0339 **swadhin: HADRONIC TAU DECAY
0340 *Here particle that are not decayed are called (except neutrinos) and their pT are added. The sum is called pT Miss.
0341 
0342 
0343         do 60 J=1,N
0344 
0345           if ((K(J,1).EQ.11).and.
0346      >        (K(J,2).EQ.15)) then      ! find tau and get it's line number
0347            idt=J
0348            idtd=K(J,4)                  ! tau decay's to what line number
0349           endif
0350 
0351           do 45 l=1,N              
0352           if  (l.EQ.idt) then      
0353                                     
0354              ptt=PYP(l,10)
0355              pxt=P(l,1)
0356              pyt=P(l,2)
0357              pht=PYP(l,16)
0358              ppt=PYP(l,14)
0359            endif
0360 45        enddo
0361 
0362           if ((K(J,1).EQ.1).and.
0363      >        (K(J,2).EQ.16)
0364      >    .and.(K(J,3).EQ.idt)) then      ! find tau neutrino and get it's line number
0365            idtnu=J
0366           endif
0367 
0368           if ((J.EQ.idtd)      ! line number is the decay of tau
0369      >   .and.(K(J,2).NE.-12)  ! this decay is not nu_ebar
0370      >   .and.(K(J,2).NE.-14)) then ! this decay is not nu_mubar
0371 
0372           ptid=0.d0
0373           pxid=0.d0
0374           pyid=0.d0
0375           pzid=0.d0
0376 
0377           do 50 I=1,N
0378 
0379           if ((K(I,1).LT.11)    ! Anything that doesn't decay = Final Product
0380      >   .and.(K(I,2).NE.12)  ! Final Product NOT AN Electron neutrino
0381      >   .and.(K(I,2).NE.14)  ! Final Product NOT AN Muon neutrino
0382      >   .and.(K(I,2).NE.16)  ! Final Product NOT AN Tau neutrino
0383      >   .and.(K(I,2).NE.-12) ! Final Product NOT AN Electron neutrino
0384      >   .and.(K(I,2).NE.-14) ! Final Product NOT AN  Muon neutrino
0385      >   .and.(K(I,2).NE.-16)) then ! Final Product NOT A Tau neutrino
0386 
0387           ptid=ptid+PYP(I,10)
0388           pxid=pxid+P(I,1)
0389           pyid=pyid+P(I,2)
0390           pzid=pzid+P(I,3)
0391           endif
0392 50        enddo
0393 
0394         !!Delta Phi < 20 cut
0395  
0396         if (pyid.GT.0.) then
0397         phimiss=acos(pxid/sqrt(pxid*pxid+pyid*pyid))*180/(3.14159)
0398         endif
0399 
0400         if (pyid.LT.0.) then
0401         phimiss=-acos(pxid/sqrt(pxid*pxid+pyid*pyid))*180/(3.14159)
0402 
0403         endif
0404 
0405         if(abs(phimiss-pht).GT.180.) then
0406         dphimiss= 360-abs(phimiss-pht)
0407         endif
0408 
0409         if(abs(phimiss-pht).GT.180.) then
0410         dphimiss= abs(phimiss-pht)
0411         endif
0412 
0413 C        write(*,*) dphimiss
0414 
0415         if (dphimiss.LT.20.) then  
0416           write(8,10)LQGPAR1(4),pxid,pyid,pxt,
0417      >    pyt,pht,ppt  
0418         endif
0419 
0420 10      FORMAT(I8,6(1PE14.6))
0421            endif
0422 60      enddo
0423 
0424 C...List first few events.
0425        LQGPAR1(4)=LQGPAR1(4)+1
0426        if(LQGPAR1(4).LE.LQGPAR1(2))         CALL PYLIST(2)
0427 C.. Swadhin
0428        if(mod(LQGPAR1(4),1000).eq.0)then
0429         write(6,1000) LQGPAR1(4)
0430 1000    format('>>>>>> ',I8,
0431      >         ' events succesfully generated    <<<<<<')
0432        endif
0433 
0434 C Write to output file
0435 C Josh
0436 
0437 C      WRITE(29,*) 0, LQGPAR1(4), 1, MSTI(1), MSTI(12), MSTI(16), 
0438 C     &           PARI(34), MSTI(15), PARI(33), PARI(53)
0439        WRITE(29,*) '   0          1    1   95 2212         21  0.000000
0440      & 21 0.000000     0.000000      0.03587194790        
0441      & 0.00016131040 0.00000022484      718.31915283203      
0442      & 382.32009887695
0443      &   0.00773132309        0.00000000000           0.000000000
0444      &   0.000000000           0.000000000           0.000000000
0445      &   0.000000000           0.000000000           0.000000000
0446      &   0.000000000           0.000000000           0.000000000
0447      &   1.422321899             670'
0448       WRITE(29,*)'============================================'
0449        DO 100 L=1,N
0450 C         WRITE(29,40) L, K(L,1), K(L,2), K(L,3)
0451 C40       FORMAT(I5.2, I5.2, I7.2, I5.2)
0452           if(L == 3) then
0453           WRITE(29,*)'============================================'
0454           endif
0455 C          if((K(L,1) .NE. 21) .AND. (KSprev .EQ. 21))
0456 C     &           WRITE(29,*)'============================================'
0457 C          WRITE(29,40) L, K(L,1), K(L,2), K(L,3), P(L,1), P(L,2)
0458 C     &                 , P(L,3), P(L,4), P(L,5)
0459 C40        FORMAT(I5.2, I5.2, I7.2, I5.2, F10.5, F10.5, F10.5
0460 C     &           , F10.5, F10.5)
0461           WRITE(29,40) L, K(L,1), K(L,2), K(L,3), K(L,4), K(L,5),
0462      &                 P(L,1), P(L,2)
0463      &                 , P(L,3), P(L,4), P(L,5), V(L,1), V(L,2), V(L,3)
0464 40        FORMAT(I5, I5, I7, I5, I5, I5, F10.5, F10.5, F10.5
0465      &           , F10.5, F10.5, F15.5, F15.5, F15.5)
0466           KSprev = K(L,1)
0467 100    continue
0468        WRITE(29,*)'=============== Event finished ==============='
0469        CALL PYHIST
0470 
0471        if(LQGPAR1(3).ne.0)then
0472 C...Fill histos
0473         CALL HF1(1000,sngl(x),1.) 
0474         if(genproc.eq.1)CALL HF1(1001,sngl(x),1.) 
0475         if(genproc.eq.2)CALL HF1(1002,sngl(x),1.) 
0476         CALL HF1(2000,sngl(y),1.) 
0477         if(genproc.eq.1)CALL HF1(2001,sngl(y),1.) 
0478         if(genproc.eq.2)CALL HF1(2002,sngl(y),1.) 
0479         CALL HF1(3000,log10(sngl(Q2)),1.) 
0480 
0481 C... final energy and momentum checks
0482         px_sum=0.
0483         py_sum=0.
0484         pz_sum=0.
0485         e_sum=0.
0486         cme=0.
0487         do 222 i=1,N
0488          if(K(I,1).le.10)then
0489           px_sum=px_sum+P(I,1)
0490           py_sum=py_sum+P(I,2)
0491           pz_sum=pz_sum+P(I,3)
0492           e_sum=e_sum+P(I,4)
0493          endif
0494 222     enddo
0495         cme=sqrt(e_sum**2-px_sum**2-py_sum**2-pz_sum**2) 
0496         call hf1(5001,sngl(px_sum),1.)
0497         call hf1(5002,sngl(py_sum),1.)
0498         call hf1(5003,sngl(pz_sum),1.)
0499         call hf1(5004,sngl(e_sum),1.)
0500         call hf1(5000,sngl(cme),1.)    
0501        endif
0502       endif
0503 
0504 C-----------------------------------------------------------------
0505 C...FLAG=2 -> Third section: Termination  
0506       if(flag.eq.2)Then
0507         write(6,*)
0508 
0509 C...Pythia final table.
0510         CALL PYSTAT(1)
0511 
0512         write(6,*)
0513 C...LQGENEP final statistics.
0514         CALL LQGPRI2
0515 
0516 C...Closing Histograms.
0517         if(LQGPAR1(3).ne.0)then
0518          Call HCDIR('//lqgenep',' ')
0519          CALL HROUT(0,ICYCLE,' ')
0520          CALL HREND('lqgenep')
0521         endif
0522       endif
0523       END
0524 *
0525 C*********************************************************************
0526  
0527       SUBROUTINE PYUPEV(ISUB,SIGEV)
0528 C-------------------------------------------
0529 C...Pythia routine for user external process
0530 C------------------------------------------- 
0531 C...Double precision and integer declarations.
0532       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0533       IMPLICIT INTEGER(I-N)
0534       INTEGER PYK,PYCHGE,PYCOMP
0535 
0536 C# LQGpar1.inc #
0537       integer echar
0538       double precision ebeam,pbeam
0539       common /LQGbeam/ebeam,pbeam,echar
0540 
0541 C# LQGpdfC.inc #
0542       double precision pdfsf(20)
0543       common /LQGpdfC/ pdfsf 
0544       
0545 C# LQGKinC.inc #
0546       double precision xmax,xmin,ymax,ymin,zmax,zmin,Q2min
0547       common /LQGKinC/ xmax,xmin,ymax,ymin,zmax,zmin,Q2min
0548 
0549 C# LQGproc.inc #
0550       double precision Mlq,G1,G2
0551       Integer LQtype,l_o,q_i,q_j
0552       common /LQGproc/ Mlq,G1,G2,LQtype,l_o,q_i,q_j 
0553 
0554 C# LQGKinV.inc # 
0555       double precision S,Srad,x,y,z,Q2
0556       common /LQGKinV/ S,Srad,x,y,z,Q2
0557 
0558 C# LQGout.inc #
0559       double precision DXSEC(3),pvalence
0560       integer q_o,q_s,genproc,genprtyp,sch,uch,gproc(8)
0561       common /LQGout/ DXSEC,pvalence,q_o,q_s,genproc,genprtyp,
0562      >sch,uch,gproc
0563 C...
0564       CHARACTER CHAF*16
0565       COMMON /PYDAT4/CHAF(500,2)  
0566       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0567       COMMON/PYUPPR/NUP,KUP(20,7),NFUP,IFUP(10,2),PUP(20,5),Q2UP(0:10)
0568       SAVE /PYDAT1/,/PYUPPR/
0569 C...Local arrays and parameters.
0570       DIMENSION XPPs(-25:25),XPPu(-25:25),XPE(-25:25),TERM(20)
0571       DATA PI/3.141592653589793D0/
0572 * conversion from pb to mb
0573       DATA CONV/1.D-9/
0574 c      DATA CONV/1.D0/
0575 
0576 C...LQGENEP parameters
0577       double precision BEAMPAR,LQGPAR3
0578       integer LQGPAR1,LQGPAR2
0579       CHARACTER*256 out_file
0580       COMMON/LQGDATA/BEAMPAR(3),LQGPAR1(10),LQGPAR2(10),LQGPAR3(20),
0581      > out_file
0582 
0583 C...LQ names according to Aachen convention
0584       Character*7 LQCHA(14)
0585       DATA LQCHA /'S_0L','S_0R','~S_0R','S_1L',
0586      >            'V_1/2L','V_1/2R','~V_1/2L',
0587      >            'V_0L','V_0R','~V_0R','V_1L',
0588      >            'S_1/2L','S_1/2R','~S_1/2L'/                  
0589 C...
0590       DATA KLQ /39/
0591 *      
0592       sigev=0.d0
0593       irej=0
0594 
0595 * sigma
0596       X=pyr(0)*(Xmax-Xmin)+Xmin
0597       Y=pyr(0)*(Ymax-Ymin)+Ymin
0598       Z=1
0599       Srad=S*z
0600       Q2=S*X*Y*Z
0601 *      
0602 C... Evaluate double differential cross section
0603       call LQGDDXS
0604 *
0605       dxdydz=(Xmax-Xmin)*(Ymax-Ymin)*(Zmax-Zmin)
0606       Sigev=(DXSEC(2)+DXSEC(3))*conv*dxdydz
0607  
0608 * fill Pythia variables for the generated process
0609 * e beam         
0610       ECM2XZ=S*X*Z
0611       ECMXZ=sqrt(ECM2XZ)
0612       NUP=5
0613       KUP(1,1)=1
0614       KUP(1,2)=(echar)*-11
0615       KUP(1,3)=0
0616       KUP(1,4)=0
0617       KUP(1,5)=0
0618       KUP(1,6)=0
0619       KUP(1,7)=0
0620       PUP(1,1)=0.
0621       PUP(1,2)=0.
0622       PUP(1,4)=Z*sqrt(S)/2.d0
0623       PUP(1,3)=PUP(1,4)            
0624       PUP(1,5)=0.            
0625 * p beam
0626       KUP(2,1)=1
0627       KUP(2,2)=q_s
0628       KUP(2,3)=0
0629       KUP(2,4)=0
0630       KUP(2,5)=0
0631       KUP(2,6)=0
0632       KUP(2,7)=0
0633       if(q_s.gt.0)then
0634        KUP(2,6)=5
0635       else
0636        KUP(2,7)=5
0637       endif
0638       PUP(2,1)=0.
0639       PUP(2,2)=0.
0640       PUP(2,4)=X*sqrt(S)/2.d0
0641       PUP(2,3)=-PUP(2,4)            
0642       PUP(2,5)=0.            
0643 * LQ
0644       KUP(3,1)=2
0645       KUP(3,2)=KLQ      
0646       CHAF(pycomp(KLQ),1)=LQCHA(LQTYPE)      
0647       KUP(3,3)=0
0648       KUP(3,4)=0
0649       KUP(3,5)=0
0650       KUP(3,6)=0
0651       KUP(3,7)=0
0652       PUP(3,1)=PUP(2,1)+PUP(1,1)
0653       PUP(3,2)=PUP(2,2)+PUP(1,2)
0654       PUP(3,3)=PUP(2,3)+PUP(1,3)            
0655       PUP(3,4)=sqrt(ECM2XZ+PUP(3,1)**2+PUP(3,2)**2+PUP(3,3)**2)            
0656       PUP(3,5)=sqrt(ECM2XZ)            
0657 
0658 * final state in sub-system cm.
0659 *   final state lepton
0660       theta=acos(1.d0-2.d0*Y)
0661       PHI=2D0*PI*PYR(0)
0662       rtshat=ECMXZ
0663       KUP(4,1)=1
0664       KUP(4,2)=echar*-(11+2*(l_o-1))
0665       KUP(5,1)=1
0666       KUP(5,2)=q_o
0667       PUP(4,5)=PYMASS(KUP(4,2))
0668       PUP(5,5)=PYMASS(KUP(5,2))
0669       PUP44=0.5D0*(RTSHAT**2+PUP(4,5)**2-PUP(5,5)**2)/RTSHAT      
0670       PUP54=RTSHAT-PUP44
0671       KUP(4,3)=3
0672       KUP(4,4)=0
0673       KUP(4,5)=0
0674       KUP(4,6)=0
0675       KUP(4,7)=0
0676       if(irej.eq.1.and.PUP44**2-PUP(4,5)**2.lt.0)then
0677        PMOD=1.d0
0678       else
0679        PMOD=sqrt(PUP44**2-PUP(4,5)**2)      
0680       endif
0681       PUP(4,1)=PMOD*sin(theta)*cos(phi)
0682       PUP(4,2)=PMOD*sin(theta)*sin(phi)
0683       PUP43=PMOD*cos(theta)      
0684       PUP44=PUP(3,5)/2.d0
0685 *   final state quark
0686       KUP(5,3)=3
0687       KUP(5,4)=0
0688       KUP(5,5)=0
0689       KUP(5,6)=0
0690       KUP(5,7)=0
0691       if(q_o.gt.0)then
0692        KUP(5,4)=2
0693       else
0694        KUP(5,5)=2
0695       endif
0696       PUP(5,1)=0.
0697       PUP(5,2)=0.
0698       if(irej.eq.1.and.PUP54**2-PUP(5,5)**2.lt.0)then
0699        PMOD=1.d0
0700       else
0701        PMOD=sqrt(PUP54**2-PUP(5,5)**2)      
0702       endif
0703       PUP(5,1)=-PUP(4,1)
0704       PUP(5,2)=-PUP(4,2)
0705       PUP53=-PUP43      
0706 *   Longitudinal boost to cm frame
0707       beta=(z-x)/(z+x)
0708       gamma=0.5d0*(z+x)/sqrt(x*z)
0709       PUP(4,3)=GAMMA*(PUP43+BETA*PUP44)
0710       PUP(4,4)=GAMMA*(PUP44+BETA*PUP43)
0711       PUP(5,3)=GAMMA*(PUP53+BETA*PUP54)
0712       PUP(5,4)=GAMMA*(PUP54+BETA*PUP53)
0713 * 
0714       NFUP=1
0715       IFUP(1,1)=4
0716       IFUP(1,2)=5
0717       Q2UP(0)=Q2
0718       Q2UP(1)=Q2
0719       RETURN
0720       END
0721 *
0722 
0723       SUBROUTINE LQGDDXS
0724 C----------------------------------------------
0725 C...Evaluate double differential cross section 
0726 C...           d^2 sigma / dx dy
0727 C----------------------------------------------
0728 *
0729       implicit none
0730 *
0731 C# LQGpar1.inc #
0732       integer echar
0733       double precision ebeam,pbeam
0734       common /LQGbeam/ebeam,pbeam,echar
0735             
0736 C# LQGpdfC.inc #
0737       double precision pdfsf(20)
0738       common /LQGpdfC/ pdfsf
0739             
0740 C# LQGKinC.inc #
0741       double precision xmax,xmin,ymax,ymin,zmax,zmin,Q2min
0742       common /LQGKinC/ xmax,xmin,ymax,ymin,zmax,zmin,Q2min
0743             
0744 C# LQGproc.inc #
0745       double precision Mlq,G1,G2
0746       Integer LQtype,l_o,q_i,q_j
0747       common /LQGproc/ Mlq,G1,G2,LQtype,l_o,q_i,q_j 
0748 
0749 C# LQGKinV.inc #
0750       double precision S,Srad,x,y,z,Q2
0751       common /LQGKinV/ S,Srad,x,y,z,Q2
0752             
0753 C# LQGout.inc #
0754       double precision DXSEC(3),pvalence
0755       integer q_o,q_s,genproc,genprtyp,sch,uch,gproc(8)
0756       common /LQGout/ DXSEC,pvalence,q_o,q_s,genproc,genprtyp,
0757      >sch,uch,gproc
0758 
0759       double precision DSIGMADXDY(4)
0760       double precision rand(1),sfrac1,sfrac2,ufrac1,ufrac2
0761       double precision pvalences_u,pvalences_d,pvalenceu_u,pvalenceu_d
0762 *
0763 cc--------------------------------------------------------
0764 C...Leptoquark types ranges from 1 to 14.
0765 C 
0766 C     1->S_0 LEFT
0767 C     2->S_0 RIGHT
0768 C     3->~S_0 RIGHT
0769 C     4->S_1 LEFT
0770 C     5->V_1/2 LEFT
0771 C     6->V_1/2 RIGHT
0772 C     7->~V_1/2 LEFT
0773 C     8->V_0 LEFT
0774 C     9->V_0 RIGHT
0775 C     10->~V_0 RIGHT
0776 C     11->V_1 LEFT
0777 C     12->S_1/2 LEFT
0778 C     13->S_1/2 RIGHT   
0779 C     14->~S_1/2 LEFT
0780 C
0781 C       DSIGMADXDY(4) - Double differential cross section
0782 C                       DSIGMADXDY(1) = Standard Model term (SM processes)
0783 C                       DSIGMADXDY(2) = Interference term between SM and LQ
0784 C                       DSIGMADXDY(3) = LQ term - u channel
0785 C                       DSIGMADXDY(4) = LQ term - s channel
0786 C ========================================================================
0787 C INPUT PARAMETERS:
0788 C                   X - standard DIS x variable 
0789 C                   Y - standard DIS y variable
0790 C                   S - Center of mass energy
0791 C                 MLQ - Leptoquark mass
0792 C                  G1 - initial state coupling
0793 C                  G2 - final state coupling
0794 C                 l_o - generation of the outcoming lepton
0795 C               echar - charge of the incoming lepton
0796 C                 q_i - generation of initial state quark
0797 C                 q_j - generation of the final state quark
0798 C              LQTYPE - Leptoquark type (see table above)
0799 C
0800 C OUTPUT PARAMETERS:
0801 C                       DXSEC = double differential cross section (pb):
0802 C                            DXSEC(1)= LQ-SM interference term
0803 C                            DXSEC(2)= LQ term - u channel 
0804 C                            DXSEC(3)= LQ term - s channel
0805 C                       q_o = output quark (from LQ decay):
0806 C                            1 down        -1 antidown 
0807 C                            2 up          -2 antiup
0808 C                            3 strange     -3 antistrange 
0809 C                            4 charm       -4 anticharm
0810 C                            5 bottom      -5 antibottom
0811 C                            6 top         -6 antitop
0812 C            
0813 C----------------------------------------------------------
0814       double precision pyr
0815       double precision C_R_P,C_L_P,C_R_E,C_L_E
0816      &                ,C_R_U,C_L_U,C_R_D,C_L_D
0817      &                ,B_RR_U,B_RL_U
0818      &                ,B_LR_U,B_LL_U
0819      &                ,B_RR_D,B_RL_D
0820      &                ,B_LR_D,B_LL_D
0821       double precision CCCf2u,DDDf2u,EEEf2u,FFFf2u,
0822      &                 CCCf2d,DDDf2d,EEEf2d,FFFf2d
0823       double precision CCCf0u,DDDf0u,EEEf0u,FFFf0u,
0824      &                 CCCf0d,DDDf0d,EEEf0d,FFFf0d
0825       double precision CCCv2u,DDDv2u,EEEv2u,FFFv2u,
0826      &                 CCCv2d,DDDv2d,EEEv2d,FFFv2d
0827       double precision CCCv0u,DDDv0u,EEEv0u,FFFv0u,
0828      &                 CCCv0d,DDDv0d,EEEv0d,FFFv0d
0829 
0830       double precision weakmix,Mz2,Gz2,Gf,alpha,A,P,pi
0831       parameter (weakmix=0.2315, Mz2=(91.187)**2,Gz2=(2.490)**2)
0832       parameter (Gf=0.0000116639,alpha=1.0/137.036, A=27.5, P=820.0)
0833       parameter (pi=3.141592653589793d0)
0834 
0835       double precision UPV,DNV,USEA,DSEA,STR,CHM,BOT,TOP,GL
0836       double precision UPQs(3),DNQs(3),UPQBs(3),DNQBs(3)
0837       double precision UPQu(3),DNQu(3),UPQBu(3),DNQBu(3)
0838       double precision scales,scaleu,LambdaL2,LambdaR2,Ms2
0839      &                ,aaa,bbb,ggg
0840      &                ,LambdaL2_1,LambdaL2_2
0841      &                ,LambdaR2_1,LambdaR2_2
0842      &                ,GAM
0843 
0844 c......LH 0, RH 1
0845       integer LRindex(14)
0846       data LRindex/0,1,1,0,0,1,0,0,1,1,0,0,1,0/
0847       integer SVindex(14)
0848       data SVindex/1,1,1,1,2,2,2,3,3,3,3,4,4,4/
0849 
0850 * protection
0851       if(y.eq.1)y=1.d0-1.d-13
0852 *
0853       DSIGMADXDY(1)=0.d0
0854       DSIGMADXDY(2)=0.d0
0855       DSIGMADXDY(3)=0.d0
0856       DSIGMADXDY(4)=0.d0
0857       q_o=0
0858       pvalence=0.
0859       pvalences_u=0.
0860       pvalenceu_u=0.
0861       pvalences_d=0.
0862       pvalenceu_d=0.
0863 *
0864       C_R_P = weakmix
0865       C_L_P = -0.5+weakmix
0866       C_R_U = -2.0*weakmix/3.0
0867       C_L_U = 0.5-2.0*weakmix/3.0
0868       C_R_D = weakmix/3.0
0869       C_L_D = -0.5+weakmix/3.0
0870       C_R_E = weakmix
0871       C_L_E = -0.5+weakmix
0872 
0873       Ms2=(MLQ)**2
0874       Q2=Srad*X*Y
0875       if(Q2.lt.Q2min)goto 999
0876 * u channel densities
0877       scaleu=sqrt(Srad*X*(1.-Y))    
0878       CALL STRUCTM(X,scaleu,UPV
0879      &               ,DNV,USEA,DSEA,STR,CHM,BOT,TOP,GL)                   
0880       if(UPV+USEA.gt.0)
0881      & pvalenceu_u=UPV/(UPV+USEA)
0882       if(DNV+DSEA.gt.0)
0883      & pvalenceu_d=DNV/(DNV+DSEA)
0884       if(echar.eq.1)then
0885 * case e+, mu+, tau+ 
0886        UPQu(1)=UPV+USEA
0887        UPQBu(1)=USEA
0888        UPQu(2)=CHM
0889        UPQBu(2)=CHM
0890        UPQu(3)=TOP
0891        UPQBu(3)=TOP
0892        DNQu(1)=DNV+DSEA
0893        DNQBu(1)=DSEA
0894        DNQu(2)=STR
0895        DNQBu(2)=STR
0896        DNQu(3)=BOT
0897        DNQBu(3)=BOT
0898       elseif(echar.eq.-1)then
0899 * case e-, mu-, tau- 
0900        UPQu(1)=USEA
0901        UPQBu(1)=UPV+USEA
0902        UPQu(2)=CHM
0903        UPQBu(2)=CHM
0904        UPQu(3)=TOP
0905        UPQBu(3)=TOP
0906        DNQu(1)=DSEA
0907        DNQBu(1)=DNV+DSEA
0908        DNQu(2)=STR
0909        DNQBu(2)=STR
0910        DNQu(3)=BOT
0911        DNQBu(3)=BOT
0912       endif          
0913 * s channel densities
0914       scales=sqrt(Srad*X)    
0915       CALL STRUCTM(X,scales,UPV
0916      &               ,DNV,USEA,DSEA,STR,CHM,BOT,TOP,GL)        
0917       if(UPV+USEA.gt.0)
0918      & pvalences_u=UPV/(UPV+USEA)
0919       if(DNV+DSEA.gt.0)
0920      & pvalences_d=DNV/(DNV+DSEA)           
0921       if(echar.eq.1)then
0922 * case e+, mu+, tau+ 
0923        UPQs(1)=UPV+USEA
0924        UPQBs(1)=USEA
0925        UPQs(2)=CHM
0926        UPQBs(2)=CHM
0927        UPQs(3)=TOP
0928        UPQBs(3)=TOP
0929        DNQs(1)=DNV+DSEA
0930        DNQBs(1)=DSEA
0931        DNQs(2)=STR
0932        DNQBs(2)=STR
0933        DNQs(3)=BOT
0934        DNQBs(3)=BOT
0935       elseif(echar.eq.-1)then
0936 * case e-, mu-, tau- 
0937        UPQs(1)=USEA
0938        UPQBs(1)=UPV+USEA
0939        UPQs(2)=CHM
0940        UPQBs(2)=CHM
0941        UPQs(3)=TOP
0942        UPQBs(3)=TOP
0943        DNQs(1)=DSEA
0944        DNQBs(1)=DNV+DSEA
0945        DNQs(2)=STR
0946        DNQBs(2)=STR
0947        DNQs(3)=BOT
0948        DNQBs(3)=BOT
0949       endif
0950 *         
0951       aaa = Q2*(Q2+Mz2)/((Q2+Mz2)**2+Mz2*Gz2)
0952       bbb = sqrt(2.0)*Gf*Mz2/(pi*alpha)
0953     
0954       B_RR_U = -2./3. + aaa*bbb*C_R_P*C_R_U
0955       B_RL_U = -2./3. + aaa*bbb*C_R_P*C_L_U
0956       B_LR_U = -2./3. + aaa*bbb*C_L_P*C_R_U
0957       B_LL_U = -2./3. + aaa*bbb*C_L_P*C_L_U
0958       B_RR_D = 1./3. + aaa*bbb*C_R_P*C_R_D
0959       B_RL_D = 1./3. + aaa*bbb*C_R_P*C_L_D
0960       B_LR_D = 1./3. + aaa*bbb*C_L_P*C_R_D
0961       B_LL_D = 1./3. + aaa*bbb*C_L_P*C_L_D
0962 
0963       IF (SVindex(LQTYPE).EQ.1) THEN
0964 
0965        CCCf2u=Q2*(1.-Y)**2*(UPV+USEA)
0966      &                /(4.0*pi*alpha)
0967        DDDf2u=Q2*USEA/(4.0*pi*alpha)
0968        EEEf2u=(Q2-X*Srad)**2*Y**2*(UPQu(q_j))
0969      &            /(64.0*(pi*alpha)**2)
0970        FFFf2u=Q2**2*UPQBs(q_i)/(64.0*(pi*alpha)**2)
0971 
0972        CCCf2d=Q2*(1.-Y)**2*(DNV+DSEA)
0973      &                /(4.0*pi*alpha)
0974        DDDf2d=Q2*DSEA/(4.0*pi*alpha)
0975        EEEf2d=(Q2-X*Srad)**2*Y**2*(DNQu(q_j))
0976      &            /(64.0*(pi*alpha)**2)
0977        FFFf2d=Q2**2*DNQBs(q_i)/(64.0*(pi*alpha)**2)
0978 
0979       ELSE IF (SVindex(LQTYPE).EQ.4) THEN
0980 
0981        CCCf0u=Q2*(1.-Y)**2*USEA/(4.0*pi*alpha)
0982        DDDf0u=Q2*(UPV+USEA)/(4.0*pi*alpha)
0983        EEEf0u=(Q2-X*Srad)**2*Y**2*UPQBu(q_j)
0984      &            /(64.0*(pi*alpha)**2)
0985        FFFf0u=Q2**2*(UPQs(q_i))/(64.0*(pi*alpha)**2)
0986 
0987        CCCf0d=Q2*(1.-Y)**2*DSEA/(4.0*pi*alpha)
0988        DDDf0d=Q2*(DNV+DSEA)/(4.0*pi*alpha)
0989        EEEf0d=(Q2-X*Srad)**2*Y**2*DNQBu(q_j)
0990      &            /(64.0*(pi*alpha)**2)
0991        FFFf0d=Q2**2*(DNQs(q_i))/(64.0*(pi*alpha)**2)
0992 
0993       ELSE IF (SVindex(LQTYPE).EQ.2) THEN
0994 
0995        CCCv2u=Q2*(1.-Y)**2*USEA/(2.0*pi*alpha)
0996        DDDv2u=Q2*(UPV+USEA)/(2.0*pi*alpha)
0997        EEEv2u=Q2**2*UPQBs(q_i)/(16.0*(pi*alpha)**2)
0998        FFFv2u=(Q2-X*Srad)**2*Y**2*(UPQu(q_j))
0999      &            /(16.0*(pi*alpha)**2*(1.0-Y)**2)
1000 
1001        CCCv2d=Q2*(1.-Y)**2*DSEA/(2.0*pi*alpha)
1002        DDDv2d=Q2*(DNV+DSEA)/(2.0*pi*alpha)
1003        EEEv2d=Q2**2*DNQBs(q_i)/(16.0*(pi*alpha)**2)
1004        FFFv2d=(Q2-X*Srad)**2*Y**2*(DNQu(q_j))
1005      &            /(16.0*(pi*alpha)**2*(1.0-Y)**2)
1006 
1007       ELSE IF (SVindex(LQTYPE).EQ.3) THEN
1008 
1009        CCCv0u=Q2*(1.-Y)**2*(UPV+USEA)
1010      &            /(2.0*pi*alpha)
1011        DDDv0u=Q2*USEA/(2.0*pi*alpha)
1012        EEEv0u=Q2**2*(UPQs(q_i))/(16.0*(pi*alpha)**2)
1013        FFFv0u=(Q2-X*Srad)**2*Y**2*UPQBu(q_j)
1014      &            /(16.0*(pi*alpha)**2*(1.0-Y)**2)
1015 
1016        CCCv0d=Q2*(1.-Y)**2*(DNV+DSEA)
1017      &              /(2.0*pi*alpha)
1018        DDDv0d=Q2*DSEA/(2.0*pi*alpha)
1019        EEEv0d=Q2**2*(DNQs(q_i))/(16.0*(pi*alpha)**2)
1020        FFFv0d=(Q2-X*Srad)**2*Y**2*DNQBu(q_j)
1021      &            /(16.0*(pi*alpha)**2*(1.0-Y)**2)
1022       ENDIF
1023 
1024       DSIGMADXDY(1)=(B_RL_U**2+B_LR_U**2+
1025      &              (B_RR_U**2+B_LL_U**2)*(1.-Y)**2)
1026      &               *(UPV+USEA+CHM+TOP)+
1027      &              (B_RL_D**2+B_LR_D**2+
1028      &              (B_RR_D**2+B_LL_D**2)*(1.-Y)**2)
1029      &               *(DNV+DSEA+STR+BOT)+
1030      &              (B_RR_U**2+B_LL_U**2+
1031      &              (B_RL_U**2+B_LR_U**2)*(1.-Y)**2)
1032      &               *(USEA+CHM+TOP) +
1033      &              (B_RR_D**2+B_LL_D**2+
1034      &              (B_RL_D**2+B_LR_D**2)*(1.-Y)**2)
1035      &               *(DSEA+STR+BOT)
1036 
1037       IF (LRindex(LQTYPE).EQ.0) THEN
1038          LambdaL2_1=G1
1039          LambdaL2_2=G2
1040          if(G2.ne.0)then
1041           LambdaL2=LambdaL2_1*LambdaL2_2
1042          else
1043           LambdaL2=G1*G1
1044          endif
1045          LambdaR2_1=0.0
1046          LambdaR2_2=0.0
1047          LambdaR2=LambdaR2_1*LambdaR2_2
1048       ELSE IF (LRindex(LQTYPE).EQ.1) THEN
1049          LambdaR2_1=G1
1050          LambdaR2_2=G2
1051          if(G2.ne.0)then
1052           LambdaR2=LambdaR2_1*LambdaR2_2
1053          else
1054           LambdaR2=G1*G1
1055          endif
1056          LambdaL2_1=0.0
1057          LambdaL2_2=0.0
1058          LambdaL2=LambdaL2_1*LambdaL2_2
1059       ENDIF
1060 
1061       IF (LQTYPE.EQ.1.or.LQTYPE.EQ.2) THEN
1062          GAM=MLQ*(sqrt(2.0)*LambdaL2_1+LambdaR2_1)**2
1063          if(l_o.gt.1)GAM=GAM+
1064      &           MLQ*(sqrt(2.0)*LambdaL2_2+LambdaR2_2)**2
1065          AAA=(X*Srad-Ms2)**2+
1066      &        Ms2*GAM**2/((16.0*pi)**2)
1067          BBB=LambdaR2*B_RR_U+LambdaL2*B_LL_U
1068          DSIGMADXDY(2)=BBB*CCCf2u/(Q2-X*Srad-Ms2)+
1069      &             BBB*DDDf2u*(X*Srad-Ms2)/AAA
1070          GGG=(LambdaR2+LambdaL2)**2
1071          if(q_i.ne.3)then
1072           DSIGMADXDY(3)=EEEf2u*GGG/(Q2-X*Srad-Ms2)**2
1073          else
1074 * Top quark in the final state
1075           DSIGMADXDY(3)=0.d0
1076          endif         
1077          if(q_j.ne.3)then
1078           DSIGMADXDY(4)=FFFf2u*GGG/AAA
1079          else
1080 * Top quark in the final state
1081           DSIGMADXDY(4)=0.d0
1082          endif
1083 * output quark
1084          if(DSIGMADXDY(3)+DSIGMADXDY(4).eq.0)then
1085           goto 999
1086          endif          
1087          sfrac1=DSIGMADXDY(4)/(DSIGMADXDY(4)+DSIGMADXDY(3))
1088 c         call rm48(rand,1)
1089          rand(1)=pyr(0)
1090          if(rand(1).lt.sfrac1)then
1091 * s channel: e+ ub -> l+ ub
1092           genproc=1
1093           genprtyp=3
1094           if(q_i.eq.1)then
1095            if(echar.lt.0)pvalence=pvalences_u
1096           endif
1097           if(q_j.eq.1)q_o=-2
1098           if(q_j.eq.2)q_o=-4
1099           if(q_j.eq.3)q_o=-6
1100           if(q_i.eq.1)q_s=-2
1101           if(q_i.eq.2)q_s=-4
1102           if(q_i.eq.3)q_s=-6          
1103          else
1104 * u channel: e+ u -> l+ u
1105           genproc=2
1106           genprtyp=5
1107           if(q_j.eq.1)then
1108            if(echar.gt.0)pvalence=pvalenceu_u
1109           endif
1110           if(q_i.eq.1)q_o=2
1111           if(q_i.eq.2)q_o=4
1112           if(q_i.eq.3)q_o=6
1113           if(q_j.eq.1)q_s=2
1114           if(q_j.eq.2)q_s=4
1115           if(q_j.eq.3)q_s=6          
1116          endif
1117       ELSE IF (LQTYPE.EQ.3) THEN
1118          GAM=MLQ*(LambdaL2_1+LambdaR2_1)**2
1119          if(l_o.gt.1)GAM=GAM+
1120      &            MLQ*(LambdaL2_2+LambdaR2_2)**2
1121          AAA=(X*Srad-Ms2)**2+
1122      &        Ms2*GAM**2/((16.0*pi)**2)
1123          BBB=LambdaR2*B_RR_D+LambdaL2*B_LL_D
1124          DSIGMADXDY(2)=BBB*CCCf2d/(Q2-X*Srad-Ms2)+
1125      &             BBB*DDDf2d*(X*Srad-Ms2)/AAA
1126          GGG=(LambdaR2+LambdaL2)**2
1127          DSIGMADXDY(3)=EEEf2d*GGG/(Q2-X*Srad-Ms2)**2
1128          DSIGMADXDY(4)=FFFf2d*GGG/AAA
1129 * output quark
1130          if(DSIGMADXDY(3)+DSIGMADXDY(4).eq.0)then
1131           goto 999
1132          endif          
1133          sfrac1=DSIGMADXDY(4)/(DSIGMADXDY(4)+DSIGMADXDY(3))
1134 c         call rm48(rand,1)
1135          rand(1)=pyr(0)
1136          if(rand(1).lt.sfrac1)then
1137 * s channel: e+ db -> l+ db
1138           genproc=1
1139           genprtyp=4
1140           if(q_i.eq.1)then
1141            if(echar.lt.0)pvalence=pvalences_d
1142           endif
1143           if(q_j.eq.1)q_o=-1
1144           if(q_j.eq.2)q_o=-3
1145           if(q_j.eq.3)q_o=-5
1146           if(q_i.eq.1)q_s=-1
1147           if(q_i.eq.2)q_s=-3
1148           if(q_i.eq.3)q_s=-5          
1149          else
1150 * u channel: e+ d -> l+ d
1151           genproc=2
1152           genprtyp=6
1153           if(q_j.eq.1)then
1154            if(echar.gt.0)pvalence=pvalenceu_d
1155           endif
1156           if(q_i.eq.1)q_o=1
1157           if(q_i.eq.2)q_o=3
1158           if(q_i.eq.3)q_o=5
1159           if(q_j.eq.1)q_s=1
1160           if(q_j.eq.2)q_s=3
1161           if(q_j.eq.3)q_s=5
1162          endif
1163       ELSE IF (LQTYPE.EQ.4) THEN
1164          GAM=MLQ*(sqrt(2.0)*LambdaL2_1+LambdaR2_1)**2
1165          if(l_o.gt.1)GAM=GAM+
1166      &            MLQ*(sqrt(2.0)*LambdaL2_2+LambdaR2_2)**2
1167          AAA=(X*Srad-Ms2)**2+
1168      &        Ms2*GAM**2/((16.0*pi)**2)
1169          BBB=LambdaR2*B_RR_D+2.0*LambdaL2*B_LL_D
1170          DSIGMADXDY(2)=BBB*CCCf2d/(Q2-X*Srad-Ms2)+
1171      &             BBB*DDDf2d*(X*Srad-Ms2)/AAA
1172          GGG=(LambdaR2+2.0*LambdaL2)**2
1173          DSIGMADXDY(3)=EEEf2d*GGG/(Q2-X*Srad-Ms2)**2
1174          DSIGMADXDY(4)=FFFf2d*GGG/AAA
1175          sfrac1=DSIGMADXDY(4)
1176          ufrac1=DSIGMADXDY(3)
1177          AAA=(X*Srad-Ms2)**2+
1178      &        Ms2*GAM**2/((16.0*pi)**2)
1179          BBB=LambdaR2*B_RR_U+LambdaL2*B_LL_U
1180          DSIGMADXDY(2)=DSIGMADXDY(2)+BBB*CCCf2u/(Q2-X*Srad-Ms2)+
1181      &             BBB*DDDf2u*(X*Srad-Ms2)/AAA
1182          GGG=(LambdaR2+LambdaL2)**2
1183 * no Top quark in the final state, otherwise no u-channel contribution
1184          if(q_i.ne.3)then
1185           DSIGMADXDY(3)=DSIGMADXDY(3)+EEEf2u*GGG/(Q2-X*Srad-Ms2)**2
1186          endif
1187 * no Top quark in the final state, otherwise no s-channel contribution
1188          if(q_j.ne.3)then
1189           DSIGMADXDY(4)=DSIGMADXDY(4)+FFFf2u*GGG/AAA
1190          endif
1191 * output quark
1192          if(DSIGMADXDY(3)+DSIGMADXDY(4).eq.0)then
1193           goto 999
1194          endif          
1195          sfrac2=DSIGMADXDY(4)-sfrac1
1196          ufrac2=DSIGMADXDY(3)-ufrac1
1197          sfrac1=sfrac1/(DSIGMADXDY(3)+DSIGMADXDY(4))
1198          sfrac2=sfrac2/(DSIGMADXDY(3)+DSIGMADXDY(4))
1199          ufrac1=ufrac1/(DSIGMADXDY(3)+DSIGMADXDY(4))
1200          ufrac2=ufrac2/(DSIGMADXDY(3)+DSIGMADXDY(4))
1201 c         call rm48(rand,1)
1202          rand(1)=pyr(0)
1203          if(rand(1).lt.sfrac1)then
1204 * s channel: e+ db -> l+ db
1205           genproc=1
1206           genprtyp=4
1207           if(q_i.eq.1)then
1208            if(echar.lt.0)pvalence=pvalences_d
1209           endif
1210           if(q_j.eq.1)q_o=-1
1211           if(q_j.eq.2)q_o=-3
1212           if(q_j.eq.3)q_o=-5
1213           if(q_i.eq.1)q_s=-1
1214           if(q_i.eq.2)q_s=-3
1215           if(q_i.eq.3)q_s=-5          
1216          elseif(rand(1).lt.(sfrac1+sfrac2))then
1217 * s channel: e+ ub -> l+ ub
1218           genproc=1
1219           genprtyp=3
1220           if(q_i.eq.1)then
1221            if(echar.lt.0)pvalence=pvalences_u
1222           endif
1223           if(q_j.eq.1)q_o=-2
1224           if(q_j.eq.2)q_o=-4
1225           if(q_j.eq.3)q_o=-6
1226           if(q_i.eq.1)q_s=-2
1227           if(q_i.eq.2)q_s=-4
1228           if(q_i.eq.3)q_s=-6
1229          elseif(rand(1).lt.(sfrac1+sfrac2+ufrac1))then
1230 * u channel: e+ d -> l+ d
1231           genproc=2
1232           genprtyp=6
1233           if(q_j.eq.1)then
1234            if(echar.gt.0)pvalence=pvalenceu_d
1235           endif
1236           if(q_i.eq.1)q_o=1
1237           if(q_i.eq.2)q_o=3
1238           if(q_i.eq.3)q_o=5
1239           if(q_j.eq.1)q_s=1
1240           if(q_j.eq.2)q_s=3
1241           if(q_j.eq.3)q_s=5          
1242          else
1243 * u channel: e+ u -> l+ u
1244           genproc=2
1245           genprtyp=5
1246           if(q_j.eq.1)then
1247            if(echar.gt.0)pvalence=pvalenceu_u
1248           endif
1249           if(q_i.eq.1)q_o=2
1250           if(q_i.eq.2)q_o=4
1251           if(q_i.eq.3)q_o=6
1252           if(q_j.eq.1)q_s=2
1253           if(q_j.eq.2)q_s=4
1254           if(q_j.eq.3)q_s=6
1255          endif
1256       ELSE IF (LQTYPE.EQ.5) THEN
1257          AAA=(X*Srad-Ms2)
1258          GAM=MLQ*(LambdaL2_1+LambdaR2_1)**2
1259          if(l_o.gt.1)GAM=GAM+
1260      &            MLQ*(LambdaL2_2+LambdaR2_2)**2
1261          BBB=LambdaR2*B_RL_D+LambdaL2*B_LR_D
1262          GGG=(LambdaL2+LambdaR2)/(24.0*pi)
1263          DSIGMADXDY(2)=BBB*CCCv2d*AAA/(AAA**2+
1264      &          (Ms2*GGG)**2)+
1265      &          BBB*DDDv2d/(Q2-X*Srad-Ms2)
1266          DSIGMADXDY(4)=EEEv2d*
1267      &         ((LambdaL2**2+LambdaR2**2)*(1.-Y)**2+
1268      &         2.0*LambdaL2*LambdaR2*Y**2)/
1269      &         (AAA**2+Ms2*GAM**2/(24.0*pi)**2)
1270          DSIGMADXDY(3)=FFFv2d*(LambdaL2**2+LambdaR2**2+
1271      &          2.0*Y**2*LambdaL2*LambdaR2)/
1272      &          (Q2-X*Srad-Ms2)**2
1273 * output quark
1274          if(DSIGMADXDY(3)+DSIGMADXDY(4).eq.0)then
1275           goto 999
1276          endif          
1277          sfrac1=DSIGMADXDY(4)/(DSIGMADXDY(3)+DSIGMADXDY(4))
1278 c         call rm48(rand,1)
1279          rand(1)=pyr(0)
1280          if(rand(1).lt.sfrac1)then
1281 * s channel: e+ db -> l+ db
1282           genproc=1
1283           genprtyp=4
1284           if(q_i.eq.1)then
1285            if(echar.lt.0)pvalence=pvalences_d
1286           endif
1287           if(q_j.eq.1)q_o=-1
1288           if(q_j.eq.2)q_o=-3
1289           if(q_j.eq.3)q_o=-5
1290           if(q_i.eq.1)q_s=-1
1291           if(q_i.eq.2)q_s=-3
1292           if(q_i.eq.3)q_s=-5
1293          else
1294 * u channel: e+ d -> l+ d
1295           genproc=2
1296           genprtyp=6
1297           if(q_j.eq.1)then
1298            if(echar.gt.0)pvalence=pvalenceu_d
1299           endif
1300           if(q_i.eq.1)q_o=1
1301           if(q_i.eq.2)q_o=3
1302           if(q_i.eq.3)q_o=5
1303           if(q_j.eq.1)q_s=1
1304           if(q_j.eq.2)q_s=3
1305           if(q_j.eq.3)q_s=5          
1306          endif
1307       ELSE IF (LQTYPE.EQ.6) THEN
1308          AAA=(X*Srad-Ms2)
1309          GAM=MLQ*(LambdaL2_1+LambdaR2_1)**2
1310          if(l_o.gt.1)GAM=GAM+
1311      &            MLQ*(LambdaL2_2+LambdaR2_2)**2
1312          BBB=LambdaR2*B_RL_D+LambdaL2*B_LR_D
1313          GGG=(LambdaL2+LambdaR2)/(24.0*pi)
1314          DSIGMADXDY(2)=BBB*CCCv2d*AAA/(AAA**2+
1315      &          (Ms2*GGG)**2)+
1316      &          BBB*DDDv2d/(Q2-X*Srad-Ms2)
1317          DSIGMADXDY(4)=EEEv2d*
1318      &         ((LambdaL2**2+LambdaR2**2)*(1.-Y)**2+
1319      &         2.0*LambdaL2*LambdaR2*Y**2)/
1320      &         (AAA**2+Ms2*GAM**2/(24.0*pi)**2)
1321          DSIGMADXDY(3)=FFFv2d*(LambdaL2**2+LambdaR2**2+
1322      &          2.0*Y**2*LambdaL2*LambdaR2)/
1323      &          (Q2-X*Srad-Ms2)**2
1324 * output quark
1325          sfrac1=DSIGMADXDY(4)
1326          ufrac1=DSIGMADXDY(3)
1327          BBB=LambdaR2*B_RL_U+LambdaL2*B_LR_U
1328          GGG=(LambdaL2+LambdaR2)/(24.0*pi)
1329          DSIGMADXDY(2)=DSIGMADXDY(2)+BBB*CCCv2u*AAA/(AAA**2+
1330      &          (Ms2*GGG)**2)+
1331      &          BBB*DDDv2u/(Q2-X*Srad-Ms2)
1332 * no Top quark in the final state, otherwise no s-channel contribution
1333          if(q_j.ne.3)then
1334           DSIGMADXDY(4)=DSIGMADXDY(4)+EEEv2u*
1335      &         ((LambdaL2**2+LambdaR2**2)*(1.-Y)**2+
1336      &         2.0*LambdaL2*LambdaR2*Y**2)/
1337      &         (AAA**2+Ms2*GAM**2/(24.0*pi)**2)
1338          endif
1339 * no Top quark in the final state, otherwise no u-channel contribution
1340          if(q_i.ne.3)then
1341           DSIGMADXDY(3)=DSIGMADXDY(3)+FFFv2u*(LambdaL2**2+LambdaR2**2+
1342      &          2.0*Y**2*LambdaL2*LambdaR2)/
1343      &          (Q2-X*Srad-Ms2)**2
1344          endif
1345 * output quark
1346          if(DSIGMADXDY(3)+DSIGMADXDY(4).eq.0)then
1347           goto 999
1348          endif          
1349          sfrac2=DSIGMADXDY(4)-sfrac1
1350          ufrac2=DSIGMADXDY(3)-ufrac1
1351          sfrac1=sfrac1/(DSIGMADXDY(3)+DSIGMADXDY(4))
1352          sfrac2=sfrac2/(DSIGMADXDY(3)+DSIGMADXDY(4))
1353          ufrac1=ufrac1/(DSIGMADXDY(3)+DSIGMADXDY(4))
1354          ufrac2=ufrac2/(DSIGMADXDY(3)+DSIGMADXDY(4))
1355 c         call rm48(rand,1)
1356          rand(1)=pyr(0)
1357          if(rand(1).lt.sfrac1)then
1358 * s channel: e+ db -> l+ db
1359           genproc=1
1360           genprtyp=4
1361           if(q_i.eq.1)then
1362            if(echar.lt.0)pvalence=pvalences_d
1363           endif
1364           if(q_j.eq.1)q_o=-1
1365           if(q_j.eq.2)q_o=-3
1366           if(q_j.eq.3)q_o=-5
1367           if(q_i.eq.1)q_s=-1
1368           if(q_i.eq.2)q_s=-3
1369           if(q_i.eq.3)q_s=-5          
1370          elseif(rand(1).lt.(sfrac1+sfrac2))then
1371 * s channel: e+ ub -> l+ ub
1372           genproc=1
1373           genprtyp=3
1374           if(q_i.eq.1)then
1375            if(echar.lt.0)pvalence=pvalences_u
1376           endif
1377           if(q_j.eq.1)q_o=-2
1378           if(q_j.eq.2)q_o=-4
1379           if(q_j.eq.3)q_o=-6
1380           if(q_i.eq.1)q_s=-2
1381           if(q_i.eq.2)q_s=-4
1382           if(q_i.eq.3)q_s=-6          
1383          elseif(rand(1).lt.(sfrac1+sfrac2+ufrac1))then
1384 * u channel: e+ d -> l+ d
1385           genproc=2
1386           genprtyp=6
1387           if(q_j.eq.1)then
1388            if(echar.gt.0)pvalence=pvalenceu_d
1389           endif
1390           if(q_i.eq.1)q_o=1
1391           if(q_i.eq.2)q_o=3
1392           if(q_i.eq.3)q_o=5
1393           if(q_j.eq.1)q_s=1
1394           if(q_j.eq.2)q_s=3
1395           if(q_j.eq.3)q_s=5          
1396          else
1397 * u channel: e+ u -> l+ u
1398           genproc=2
1399           genprtyp=5
1400           if(q_j.eq.1)then
1401            if(echar.gt.0)pvalence=pvalenceu_u
1402           endif
1403           if(q_i.eq.1)q_o=2
1404           if(q_i.eq.2)q_o=4
1405           if(q_i.eq.3)q_o=6
1406           if(q_j.eq.1)q_s=2
1407           if(q_j.eq.2)q_s=4
1408           if(q_j.eq.3)q_s=6
1409          endif
1410       ELSE IF (LQTYPE.EQ.7) THEN
1411          AAA=(X*Srad-Ms2)
1412          GAM=MLQ*(LambdaL2_1+LambdaR2_1)**2
1413          if(l_o.gt.1)GAM=GAM+
1414      &            MLQ*(LambdaL2_2+LambdaR2_2)**2
1415          BBB=LambdaR2*B_RL_U+LambdaL2*B_LR_U
1416          GGG=(LambdaL2+LambdaR2)/(24.0*pi)
1417          DSIGMADXDY(2)=BBB*CCCv2u*AAA/(AAA**2+
1418      &          (Ms2*GGG)**2)+
1419      &          BBB*DDDv2u/(Q2-X*Srad-Ms2)
1420 * no Top quark in the final state, otherwise no s-channel contribution
1421          if(q_j.ne.3)then
1422           DSIGMADXDY(4)=EEEv2u*
1423      &         ((LambdaL2**2+LambdaR2**2)*(1.-Y)**2+
1424      &         2.0*LambdaL2*LambdaR2*Y**2)/
1425      &         (AAA**2+Ms2*GAM**2/(24.0*pi)**2)
1426          else    
1427           DSIGMADXDY(4)=0.d0
1428          endif
1429 * no Top quark in the final state, otherwise no s-channel contribution
1430          if(q_i.ne.3)then
1431           DSIGMADXDY(3)=FFFv2u*(LambdaL2**2+LambdaR2**2+
1432      &          2.0*Y**2*LambdaL2*LambdaR2)/
1433      &          (Q2-X*Srad-Ms2)**2
1434          else    
1435           DSIGMADXDY(3)=0.d0
1436          endif
1437 * output quark
1438          if(DSIGMADXDY(3)+DSIGMADXDY(4).eq.0)then
1439           goto 999
1440          endif          
1441          sfrac1=DSIGMADXDY(4)/(DSIGMADXDY(3)+DSIGMADXDY(4))
1442 c         call rm48(rand,1)
1443          rand(1)=pyr(0)
1444          if(rand(1).lt.sfrac1)then
1445 * s channel: e+ ub -> l+ ub
1446           genproc=1
1447           genprtyp=3
1448           if(q_i.eq.1)then
1449            if(echar.lt.0)pvalence=pvalences_u
1450           endif
1451           if(q_j.eq.1)q_o=-2
1452           if(q_j.eq.2)q_o=-4
1453           if(q_j.eq.3)q_o=-6
1454           if(q_i.eq.1)q_s=-2
1455           if(q_i.eq.2)q_s=-4
1456           if(q_i.eq.3)q_s=-6          
1457          else
1458 * u channel: e+ u -> l+ u
1459           genproc=2
1460           genprtyp=5
1461           if(q_j.eq.1)then
1462            if(echar.gt.0)pvalence=pvalenceu_u
1463           endif
1464           if(q_i.eq.1)q_o=2
1465           if(q_i.eq.2)q_o=4
1466           if(q_i.eq.3)q_o=6
1467           if(q_j.eq.1)q_s=2
1468           if(q_j.eq.2)q_s=4
1469           if(q_j.eq.3)q_s=6          
1470          endif
1471       ELSE IF (LQTYPE.EQ.8.or.LQTYPE.EQ.9) THEN
1472          AAA=(X*Srad-Ms2)
1473          GAM=MLQ*(sqrt(2.0)*LambdaL2_1+LambdaR2_1)**2
1474          if(l_o.gt.1)GAM=GAM+
1475      &            MLQ*(sqrt(2.0)*LambdaL2_2+LambdaR2_2)**2
1476          BBB=-LambdaR2*B_RR_D-LambdaL2*B_LL_D
1477          GGG=(2.0*LambdaL2+LambdaR2)/(24.0*pi)
1478          DSIGMADXDY(2)=BBB*CCCv0d*AAA/(AAA**2+
1479      &          (Ms2*GGG)**2)+
1480      &          BBB*DDDv0d/(Q2-X*Srad-Ms2)
1481          DSIGMADXDY(4)=EEEv0d*
1482      &         ((LambdaL2**2+LambdaR2**2)*(1.-Y)**2+
1483      &         2.0*LambdaL2*LambdaR2*Y**2)/
1484      &         (AAA**2+Ms2*GAM**2/(24.0*pi)**2)
1485          DSIGMADXDY(3)=FFFv0d*(LambdaL2**2+LambdaR2**2+
1486      &          2.0*Y**2*LambdaL2*LambdaR2)/
1487      &          (Q2-X*Srad-Ms2)**2
1488 * output quark
1489          if(DSIGMADXDY(3)+DSIGMADXDY(4).eq.0)then
1490           goto 999
1491          endif          
1492          sfrac1=DSIGMADXDY(4)/(DSIGMADXDY(3)+DSIGMADXDY(4))
1493 c         call rm48(rand,1)
1494          rand(1)=pyr(0)
1495          if(rand(1).lt.sfrac1)then
1496 * s channel: e+ d -> l+ d
1497           genproc=1
1498           genprtyp=2
1499           if(q_i.eq.1)then
1500            if(echar.gt.0)pvalence=pvalences_d
1501           endif
1502           if(q_j.eq.1)q_o=1
1503           if(q_j.eq.2)q_o=3
1504           if(q_j.eq.3)q_o=5
1505           if(q_i.eq.1)q_s=1
1506           if(q_i.eq.2)q_s=3
1507           if(q_i.eq.3)q_s=5
1508          else
1509 * u channel: e+ db -> l+ db
1510           genproc=2
1511           genprtyp=8
1512           if(q_j.eq.1)then
1513            if(echar.lt.0)pvalence=pvalenceu_d
1514           endif
1515           if(q_i.eq.1)q_o=-1
1516           if(q_i.eq.2)q_o=-3
1517           if(q_i.eq.3)q_o=-5
1518           if(q_j.eq.1)q_s=-1
1519           if(q_j.eq.2)q_s=-3
1520           if(q_j.eq.3)q_s=-5
1521          endif
1522       ELSE IF (LQTYPE.EQ.10) THEN
1523          AAA=(X*Srad-Ms2)
1524          GAM=MLQ*(sqrt(2.0)*LambdaL2_1+LambdaR2_1)**2
1525          if(l_o.gt.1)GAM=GAM+
1526      &            MLQ*(sqrt(2.0)*LambdaL2_2+LambdaR2_2)**2
1527          BBB=-LambdaR2*B_RR_U-2.0*LambdaL2*B_LL_U
1528          GGG=(2.0*LambdaL2+LambdaR2)/(24.0*pi)
1529          DSIGMADXDY(2)=BBB*CCCv0u*AAA/(AAA**2+
1530      &          (Ms2*GGG)**2)+
1531      &          BBB*DDDv0u/(Q2-X*Srad-Ms2)
1532          if(q_j.ne.3)then
1533           DSIGMADXDY(4)=EEEv0u*
1534      &         ((4.0*LambdaL2**2+LambdaR2**2)*(1.-Y)**2+
1535      &         2.0*2.0*LambdaL2*LambdaR2*Y**2)/
1536      &         (AAA**2+Ms2*GAM**2/(24.0*pi)**2)
1537          else
1538 * Top quark in the final state: no s-channel contribution
1539           DSIGMADXDY(4)=0.d0
1540          endif
1541          if(q_i.ne.3)then
1542           DSIGMADXDY(3)=
1543      &          FFFv0u*(4.0*LambdaL2**2+LambdaR2**2+
1544      &          2.0*Y**2*2.0*LambdaL2*LambdaR2)/
1545      &          (Q2-X*Srad-Ms2)**2
1546          else
1547 * Top quark in the final state: no u-channel contribution
1548           DSIGMADXDY(3)=0.d0
1549          endif
1550 * output quark
1551          if(DSIGMADXDY(3)+DSIGMADXDY(4).eq.0)then
1552           goto 999
1553          endif          
1554          sfrac1=DSIGMADXDY(4)/(DSIGMADXDY(3)+DSIGMADXDY(4))
1555 c         call rm48(rand,1)
1556          rand(1)=pyr(0)
1557          if(rand(1).lt.sfrac1)then
1558 * s channel: e+ u -> l+ u
1559           genproc=1
1560           genprtyp=1
1561           if(q_i.eq.1)then
1562            if(echar.gt.0)pvalence=pvalences_u
1563           endif
1564           if(q_j.eq.1)q_o=2
1565           if(q_j.eq.2)q_o=4
1566           if(q_j.eq.3)q_o=6
1567           if(q_i.eq.1)q_s=2
1568           if(q_i.eq.2)q_s=4
1569           if(q_i.eq.3)q_s=6
1570          else
1571 * u channel: e+ ub -> l+ ub
1572           genproc=2
1573           genprtyp=7
1574           if(q_j.eq.1)then
1575            if(echar.lt.0)pvalence=pvalenceu_u
1576           endif
1577           if(q_i.eq.1)q_o=-2
1578           if(q_i.eq.2)q_o=-4
1579           if(q_i.eq.3)q_o=-6
1580           if(q_j.eq.1)q_s=-2
1581           if(q_j.eq.2)q_s=-4
1582           if(q_j.eq.3)q_s=-6
1583          endif
1584       ELSE IF (LQTYPE.EQ.11) THEN
1585          AAA=(X*Srad-Ms2)
1586          GAM=MLQ*(sqrt(2.0)*LambdaL2_1+LambdaR2_1)**2
1587          if(l_o.gt.1)GAM=GAM+
1588      &            MLQ*(sqrt(2.0)*LambdaL2_2+LambdaR2_2)**2
1589          BBB=-LambdaR2*B_RR_D-LambdaL2*B_LL_D
1590          GGG=(2.0*LambdaL2+LambdaR2)/(24.0*pi)
1591          DSIGMADXDY(2)=BBB*CCCv0d*AAA/(AAA**2+
1592      &          (Ms2*GGG)**2)+
1593      &          BBB*DDDv0d/(Q2-X*Srad-Ms2)
1594          DSIGMADXDY(4)=EEEv0d*
1595      &         ((LambdaL2**2+LambdaR2**2)*(1.-Y)**2+
1596      &         2.0*LambdaL2*LambdaR2*Y**2)/
1597      &         (AAA**2+Ms2*GAM**2/(24.0*pi)**2)
1598          DSIGMADXDY(3)=FFFv0d*(LambdaL2**2+LambdaR2**2+
1599      &          2.0*Y**2*LambdaL2*LambdaR2)/
1600      &          (Q2-X*Srad-Ms2)**2
1601 * output quark
1602          sfrac1=DSIGMADXDY(4)
1603          ufrac1=DSIGMADXDY(3)
1604          BBB=-LambdaR2*B_RR_U-2.0*LambdaL2*B_LL_U
1605          GGG=(2.0*LambdaL2+LambdaR2)/(24.0*pi)
1606          DSIGMADXDY(2)=DSIGMADXDY(2)+BBB*CCCv0u*AAA/(AAA**2+
1607      &          (Ms2*GGG)**2)+
1608      &          BBB*DDDv0u/(Q2-X*Srad-Ms2)
1609          if(q_j.ne.3)then
1610 * no Top quark in the final state otherwise no s-channel contribution
1611           DSIGMADXDY(4)=DSIGMADXDY(4)+EEEv0u*
1612      &         ((4.0*LambdaL2**2+LambdaR2**2)*(1.-Y)**2+
1613      &         2.0*2.0*LambdaL2*LambdaR2*Y**2)/
1614      &         (AAA**2+Ms2*GAM**2/(24.0*pi)**2)
1615          endif
1616          if(q_i.ne.3)then
1617 * no Top quark in the final state otherwise no u-channel contribution
1618           DSIGMADXDY(3)=DSIGMADXDY(3)+
1619      &          FFFv0u*(4.0*LambdaL2**2+LambdaR2**2+
1620      &          2.0*Y**2*2.0*LambdaL2*LambdaR2)/
1621      &          (Q2-X*Srad-Ms2)**2
1622          endif
1623 * output quark
1624          if(DSIGMADXDY(3)+DSIGMADXDY(4).eq.0)then
1625           goto 999
1626          endif          
1627          sfrac2=DSIGMADXDY(4)-sfrac1
1628          ufrac2=DSIGMADXDY(3)-ufrac1
1629          sfrac1=sfrac1/(DSIGMADXDY(3)+DSIGMADXDY(4))
1630          sfrac2=sfrac2/(DSIGMADXDY(3)+DSIGMADXDY(4))
1631          ufrac1=ufrac1/(DSIGMADXDY(3)+DSIGMADXDY(4))
1632          ufrac2=ufrac2/(DSIGMADXDY(3)+DSIGMADXDY(4))
1633 c         call rm48(rand,1)
1634          rand(1)=pyr(0)
1635          if(rand(1).lt.sfrac1)then
1636 * s channel: e+ d -> l+ d
1637           genproc=1
1638           genprtyp=2
1639           if(q_i.eq.1)then
1640            if(echar.gt.0)pvalence=pvalences_d
1641           endif
1642           if(q_j.eq.1)q_o=1
1643           if(q_j.eq.2)q_o=3
1644           if(q_j.eq.3)q_o=5
1645           if(q_i.eq.1)q_s=1
1646           if(q_i.eq.2)q_s=3
1647           if(q_i.eq.3)q_s=5
1648          elseif(rand(1).lt.(sfrac1+sfrac2))then
1649 * s channel: e+ u -> l+ u
1650           genproc=1
1651           genprtyp=1
1652           if(q_i.eq.1)then
1653            if(echar.gt.0)pvalence=pvalences_u
1654           endif
1655           if(q_j.eq.1)q_o=2
1656           if(q_j.eq.2)q_o=4
1657           if(q_j.eq.3)q_o=6
1658           if(q_i.eq.1)q_s=2
1659           if(q_i.eq.2)q_s=4
1660           if(q_i.eq.3)q_s=6
1661          elseif(rand(1).lt.(sfrac1+sfrac2+ufrac1))then
1662 * u channel: e+ db -> l+ db
1663           genproc=2
1664           genprtyp=8
1665           if(q_j.eq.1)then
1666            if(echar.lt.0)pvalence=pvalenceu_d
1667           endif
1668           if(q_i.eq.1)q_o=-1
1669           if(q_i.eq.2)q_o=-3
1670           if(q_i.eq.3)q_o=-5
1671           if(q_j.eq.1)q_s=-1
1672           if(q_j.eq.2)q_s=-3
1673           if(q_j.eq.3)q_s=-5
1674          else
1675 * u channel: e+ ub -> l+ ub
1676           genproc=2
1677           genprtyp=7
1678           if(q_j.eq.1)then
1679            if(echar.lt.0)pvalence=pvalenceu_u
1680           endif
1681           if(q_i.eq.1)q_o=-2
1682           if(q_i.eq.2)q_o=-4
1683           if(q_i.eq.3)q_o=-6
1684           if(q_j.eq.1)q_s=-2
1685           if(q_j.eq.2)q_s=-4
1686           if(q_j.eq.3)q_s=-6
1687          endif
1688       ELSE IF (LQTYPE.EQ.12) THEN
1689          GAM=MLQ*(LambdaL2_1+LambdaR2_1)**2
1690          if(l_o.gt.1)GAM=GAM+
1691      &            MLQ*(LambdaL2_2+LambdaR2_2)**2
1692          AAA=(X*Srad-Ms2)**2+
1693      &        Ms2*GAM**2/((16.0*pi)**2)
1694          BBB=-LambdaR2*B_RL_U-LambdaL2*B_LR_U
1695          DSIGMADXDY(2)=BBB*CCCf0u/(Q2-X*Srad-Ms2)+
1696      &             BBB*DDDf0u*(X*Srad-Ms2)/AAA
1697          GGG=(LambdaR2+LambdaL2)**2
1698          if(q_i.ne.3)then
1699           DSIGMADXDY(3)=EEEf0u*GGG/(Q2-X*Srad-Ms2)**2
1700          else
1701 * Top quark in the final state no u-channel contribution
1702           DSIGMADXDY(3)=0.d0
1703          endif
1704          if(q_j.ne.3)then
1705           DSIGMADXDY(4)=FFFf0u*GGG/AAA
1706          else
1707 * Top quark in the final state no s-channel contribution
1708           DSIGMADXDY(4)=0.d0
1709          endif
1710 * output quark
1711          if(DSIGMADXDY(3)+DSIGMADXDY(4).eq.0)then
1712           goto 999
1713          endif          
1714          if(DSIGMADXDY(3).gt.0)then
1715           sfrac1=DSIGMADXDY(4)/(DSIGMADXDY(3)+DSIGMADXDY(4))
1716          else
1717           sfrac1=0.d0
1718          endif
1719 c         call rm48(rand,1)
1720          rand(1)=pyr(0)
1721          if(rand(1).lt.sfrac1)then
1722 * s channel: e+ u -> l+ u
1723           genproc=1
1724           genprtyp=1
1725           if(q_i.eq.1)then
1726            if(echar.gt.0)pvalence=pvalences_u
1727           endif
1728           if(q_j.eq.1)q_o=2
1729           if(q_j.eq.2)q_o=4
1730           if(q_j.eq.3)q_o=6
1731           if(q_i.eq.1)q_s=2
1732           if(q_i.eq.2)q_s=4
1733           if(q_i.eq.3)q_s=6
1734          else
1735 * u channel: e+ ub -> l+ ub
1736           genproc=2
1737           genprtyp=7
1738           if(q_j.eq.1)then
1739            if(echar.lt.0)pvalence=pvalenceu_u
1740           endif
1741           if(q_i.eq.1)q_o=-2
1742           if(q_i.eq.2)q_o=-4
1743           if(q_i.eq.3)q_o=-6
1744           if(q_j.eq.1)q_s=-2
1745           if(q_j.eq.2)q_s=-4
1746           if(q_j.eq.3)q_s=-6
1747          endif
1748       ELSE IF (LQTYPE.EQ.13) THEN
1749          GAM=MLQ*(LambdaL2_1+LambdaR2_1)**2
1750          if(l_o.gt.1)GAM=GAM+
1751      &            MLQ*(LambdaL2_2+LambdaR2_2)**2
1752          AAA=(X*Srad-Ms2)**2+
1753      &        Ms2*GAM**2/((16.0*pi)**2)
1754          BBB=-LambdaR2*B_RL_U-LambdaL2*B_LR_U
1755          DSIGMADXDY(2)=BBB*CCCf0u/(Q2-X*Srad-Ms2)+
1756      &             BBB*DDDf0u*(X*Srad-Ms2)/AAA
1757          GGG=(LambdaR2+LambdaL2)**2
1758          if(q_i.ne.3)then
1759           DSIGMADXDY(3)=EEEf0u*GGG/(Q2-X*Srad-Ms2)**2
1760          else
1761 * Top quark in the final state no u-channel contribution 
1762           DSIGMADXDY(3)=0.d0
1763          endif
1764          if(q_j.ne.3)then
1765           DSIGMADXDY(4)=FFFf0u*GGG/AAA
1766          else
1767 * Top quark in the final state no s-channel contribution 
1768           DSIGMADXDY(3)=0.d0
1769          endif
1770 * output quark
1771          sfrac1=DSIGMADXDY(4)
1772          ufrac1=DSIGMADXDY(3)
1773          BBB=-LambdaR2*B_RL_D-LambdaL2*B_LR_D
1774          DSIGMADXDY(2)=DSIGMADXDY(2)+BBB*CCCf0d/(Q2-X*Srad-Ms2)+
1775      &             BBB*DDDf0d*(X*Srad-Ms2)/AAA
1776          GGG=(LambdaR2+LambdaL2)**2
1777          DSIGMADXDY(3)=DSIGMADXDY(3)+EEEf0d*GGG/(Q2-X*Srad-Ms2)**2
1778          DSIGMADXDY(4)=DSIGMADXDY(4)+FFFf0d*GGG/AAA
1779          if(DSIGMADXDY(3)+DSIGMADXDY(4).eq.0)then
1780           goto 999
1781          endif          
1782          sfrac2=DSIGMADXDY(4)-sfrac1
1783          ufrac2=DSIGMADXDY(3)-ufrac1
1784          sfrac1=sfrac1/(DSIGMADXDY(3)+DSIGMADXDY(4))
1785          sfrac2=sfrac2/(DSIGMADXDY(3)+DSIGMADXDY(4))
1786          ufrac1=ufrac1/(DSIGMADXDY(3)+DSIGMADXDY(4))
1787          ufrac2=ufrac2/(DSIGMADXDY(3)+DSIGMADXDY(4))
1788 c         call rm48(rand,1)
1789          rand(1)=pyr(0)
1790          if(rand(1).lt.sfrac1)then
1791 * s channel: e+ u -> l+ u
1792           genproc=1
1793           genprtyp=1
1794           if(q_i.eq.1)then
1795            if(echar.gt.0)pvalence=pvalences_u
1796           endif
1797           if(q_j.eq.1)q_o=2
1798           if(q_j.eq.2)q_o=4
1799           if(q_j.eq.3)q_o=6
1800           if(q_i.eq.1)q_s=2
1801           if(q_i.eq.2)q_s=4
1802           if(q_i.eq.3)q_s=6
1803          elseif(rand(1).lt.(sfrac1+sfrac2))then
1804 * s channel: e+ d -> l+ d
1805           genproc=1
1806           genprtyp=2
1807           if(q_i.eq.1)then
1808            if(echar.gt.0)pvalence=pvalences_d
1809           endif
1810           if(q_j.eq.1)q_o=1
1811           if(q_j.eq.2)q_o=3
1812           if(q_j.eq.3)q_o=5
1813           if(q_i.eq.1)q_s=1
1814           if(q_i.eq.2)q_s=3
1815           if(q_i.eq.3)q_s=5
1816          elseif(rand(1).lt.(sfrac1+sfrac2+ufrac1))then
1817 * u channel: e+ ub -> l+ ub
1818           genproc=2
1819           genprtyp=7
1820           if(q_j.eq.1)then
1821            if(echar.lt.0)pvalence=pvalenceu_u
1822           endif
1823           if(q_i.eq.1)q_o=-2
1824           if(q_i.eq.2)q_o=-4
1825           if(q_i.eq.3)q_o=-5
1826           if(q_j.eq.1)q_s=-2
1827           if(q_j.eq.2)q_s=-4
1828           if(q_j.eq.3)q_s=-5
1829          else
1830 * u channel: e+ db -> l+ db
1831           genproc=2
1832           genprtyp=8
1833           if(q_j.eq.1)then
1834            if(echar.lt.0)pvalence=pvalenceu_d
1835           endif
1836           if(q_i.eq.1)q_o=-1
1837           if(q_i.eq.2)q_o=-3
1838           if(q_i.eq.3)q_o=-5
1839           if(q_j.eq.1)q_s=-1
1840           if(q_j.eq.2)q_s=-3
1841           if(q_j.eq.3)q_s=-5
1842          endif
1843       ELSE IF (LQTYPE.EQ.14) THEN
1844          GAM=MLQ*(LambdaL2_1+LambdaR2_1)**2
1845          if(l_o.gt.1)GAM=GAM+
1846      &            MLQ*(LambdaL2_2+LambdaR2_2)**2
1847          AAA=(X*Srad-Ms2)**2+
1848      &        Ms2*GAM**2/((16.0*pi)**2)
1849          BBB=-LambdaR2*B_RL_D-LambdaL2*B_LR_D
1850          DSIGMADXDY(2)=BBB*CCCf0d/(Q2-X*Srad-Ms2)+
1851      &             BBB*DDDf0d*(X*Srad-Ms2)/AAA
1852          GGG=(LambdaR2+LambdaL2)**2
1853          DSIGMADXDY(3)=EEEf0d*GGG/(Q2-X*Srad-Ms2)**2
1854          DSIGMADXDY(4)=FFFf0d*GGG/AAA
1855 * output quark
1856          if(DSIGMADXDY(3)+DSIGMADXDY(4).eq.0)then
1857           goto 999
1858          endif          
1859          sfrac1=DSIGMADXDY(4)/(DSIGMADXDY(3)+DSIGMADXDY(4))
1860 c         call rm48(rand,1)
1861          rand(1)=pyr(0)
1862          if(rand(1).lt.sfrac1)then
1863 * s channel: e+ d -> l+ d
1864           genproc=1
1865           genprtyp=2
1866           if(q_i.eq.1)then
1867            if(echar.gt.0)pvalence=pvalences_d
1868           endif
1869           if(q_j.eq.1)q_o=1
1870           if(q_j.eq.2)q_o=3
1871           if(q_j.eq.3)q_o=5
1872           if(q_i.eq.1)q_s=1
1873           if(q_i.eq.2)q_s=3
1874           if(q_i.eq.3)q_s=5
1875          else
1876 * u channel: e+ db -> l+ db
1877           genproc=2
1878           genprtyp=8
1879           if(q_j.eq.1)then
1880            if(echar.lt.0)pvalence=pvalenceu_d
1881           endif
1882           if(q_i.eq.1)q_o=-1
1883           if(q_i.eq.2)q_o=-3
1884           if(q_i.eq.3)q_o=-5
1885           if(q_j.eq.1)q_s=-1
1886           if(q_j.eq.2)q_s=-3
1887           if(q_j.eq.3)q_s=-5
1888          endif
1889       ENDIF
1890       DSIGMADXDY(1)=DSIGMADXDY(1)*pi*alpha**2/(Srad*(X*Y)**2)
1891       DSIGMADXDY(2)=DSIGMADXDY(2)*pi*alpha**2/(Srad*(X*Y)**2)
1892       DSIGMADXDY(3)=DSIGMADXDY(3)*pi*alpha**2/(Srad*(X*Y)**2)
1893       DSIGMADXDY(4)=DSIGMADXDY(4)*pi*alpha**2/(Srad*(X*Y)**2)
1894       DSIGMADXDY(1)=DSIGMADXDY(1)*0.38938*(10.0)**9
1895       DSIGMADXDY(2)=DSIGMADXDY(2)*0.38938*(10.0)**9
1896       DSIGMADXDY(3)=DSIGMADXDY(3)*0.38938*(10.0)**9
1897       DSIGMADXDY(4)=DSIGMADXDY(4)*0.38938*(10.0)**9
1898 *
1899 999   continue
1900 *
1901       if(echar.lt.0)q_o=-q_o
1902       if(echar.lt.0)q_s=-q_s
1903       if(l_o.eq.1)then
1904 * if lepton generation = 1 => load interference term  
1905        DXSEC(1)=DSIGMADXDY(2)
1906       else
1907        DXSEC(1)=0.d0
1908       endif
1909       DXSEC(2)=DSIGMADXDY(3)
1910       DXSEC(3)=DSIGMADXDY(4)
1911 * 
1912       RETURN
1913       END
1914 *
1915       subroutine LQGBAN
1916 C-------------------------
1917 C...Print LQGENEP banner 
1918 C-------------------------
1919       implicit none
1920       
1921       write(6,*)
1922       write(6,*) ' *************************************************'
1923       write(6,*) '-                                                 -'
1924       write(6,*) '-                      LQGENEP                    -'
1925       write(6,*) '-                                                 -'
1926       write(6,*) '-                                                 -'
1927       write(6,*) '-                                                 -'
1928       write(6,*) '-    LeptoQuark GENerator for Electron-Proton     -'
1929       write(6,*) '-                   scattering                    -'      
1930       write(6,*) '-                                                 -'
1931       write(6,*) '-                                                 -'
1932       write(6,*) '-                                                 -'
1933       write(6,*) '- Author: L.Bellagamba                            -'
1934       write(6,*) '- e-mail lorenzo.bellagamba@bo.infn.it            -'
1935       write(6,*) '-                                                 -'
1936       write(6,*) '- Version: 1.0                                    -'
1937       write(6,*) '- Date: 01.03.2001                                -'
1938       write(6,*) ' *************************************************'
1939       write(6,*)
1940       write(6,*)
1941       write(6,*)
1942       return
1943       end
1944 *
1945       subroutine LQGPRI1
1946 C------------------------
1947 C...Print Run requests
1948 C------------------------
1949 *
1950       implicit none
1951 *
1952 C...LQGENEP parameters
1953       double precision BEAMPAR,LQGPAR3
1954       integer LQGPAR1,LQGPAR2
1955       CHARACTER*256 out_file
1956       COMMON/LQGDATA/BEAMPAR(3),LQGPAR1(10),LQGPAR2(10),LQGPAR3(20),
1957      > out_file
1958 
1959 *       
1960       write(6,*) '>>>>>>>>>>>>>>  LQGENEP RUN REQUEST  <<<<<<<<<<<<<<'
1961       write(6,*)
1962       write(6,101) LQGPAR1(1)
1963       write(6,1011) BEAMPAR(1)
1964       write(6,1012) BEAMPAR(2)
1965       write(6,1013) BEAMPAR(3)
1966       write(6,102) LQGPAR2(1)
1967       write(6,103) sngl(LQGPAR3(1))
1968       write(6,104) sngl(LQGPAR3(2))
1969       write(6,105) sngl(LQGPAR3(3))
1970       write(6,106) LQGPAR2(2)
1971       write(6,107) LQGPAR2(3)
1972       write(6,108) LQGPAR2(4)
1973       write(6,*)
1974       write(6,109) sngl(LQGPAR3(4)),sngl(LQGPAR3(5))
1975       write(6,110) sngl(LQGPAR3(6)),sngl(LQGPAR3(7))
1976       write(6,111) sngl(LQGPAR3(8))
1977       write(6,*)
1978       write(6,*) '>>>>>>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<'
1979 101   format('   Number of events.................. ',I8)
1980 1011  format('   Electron beam charge.............. ',F3.0)
1981 1012  format('   Electron beam energy.............. ',F7.2)
1982 1013  format('   Proton beam energy................ ',F7.2)
1983 102   format('   LQ type...........................  ',I2)
1984 103   format('   LQ mass (GeV)..................... ',F7.2)
1985 104   format('   LQ production coupling (s-ch.).... ',F7.4)
1986 105   format('   LQ decay coupling (s-ch.)......... ',F7.4)
1987 106   format('   struck quark generation (s-ch.)... ',I2)
1988 107   format('   output quark generation (s-ch.)... ',I2)
1989 108   format('   output lepton generation.......... ',I2)
1990 109   format('   X generation range................ ',F6.3,' - ',F6.3)
1991 110   format('   Y generation range................ ',F6.3,' - ',F6.3)
1992 111   format('   minimum allowed Q2 (GeV^2)........ ',F7.2)
1993       return
1994       end
1995 *
1996       subroutine LQGPRI2
1997 C------------------------------------------
1998 C... Print LQGENEP final statistics 
1999 C------------------------------------------
2000 *
2001       Implicit none
2002 *
2003       double precision DXSEC(3),pvalence
2004       integer q_o,q_s,genproc,genprtyp,sch,uch,gproc(8)
2005       common /LQGout/ DXSEC,pvalence,q_o,q_s,genproc,genprtyp,
2006      >sch,uch,gproc
2007 *
2008 C...LQGENEP run setup parameters
2009       double precision BEAMPAR,LQGPAR3
2010       integer LQGPAR1,LQGPAR2
2011       CHARACTER*256 out_file
2012       COMMON/LQGDATA/BEAMPAR(3),LQGPAR1(10),LQGPAR2(10),LQGPAR3(20),
2013      > out_file
2014 
2015 
2016 C...LQGENEP event informations
2017       double precision LQGKPAR,LQGDDX
2018       integer LQGPID
2019       COMMON/LQGEVT/LQGKPAR(3),LQGDDX(3),LQGPID(3)
2020 
2021 C# LQGproc.inc #
2022       double precision Mlq,G1,G2
2023       Integer LQtype,l_o,q_i,q_j
2024       common /LQGproc/ Mlq,G1,G2,LQtype,l_o,q_i,q_j  
2025 
2026 C# LQGpar1.inc #
2027       integer echar
2028       double precision ebeam,pbeam
2029       common /LQGbeam/ebeam,pbeam,echar           
2030 
2031       Character*40 CHANS1,CHANS2,CHANU1,CHANU2
2032 *       
2033       CALL LQGCHAN(CHANS1,CHANS2,CHANU1,CHANU2)
2034 *
2035       write(6,*) '>>>>>>>>>>>  LQGENEP FINAL STATISTICS  <<<<<<<<<<<'
2036       write(6,*)
2037       write(6,101) LQGPAR1(4)
2038       write(6,104) sch
2039       write(6,106) uch
2040       write(6,*)
2041       write(6,*)
2042       write(6,*) '            Details of the generation             '
2043       write(6,*) '           ---------------------------            '
2044       write(6,*)
2045       write(6,*) '=================================================='
2046       write(6,*) 'I                                       I        I'
2047       write(6,*) 'I                process                I events I'
2048       write(6,*) 'I                                       I        I'
2049       write(6,*) '=================================================='
2050       write(6,*) 'I                                       I        I'      
2051       if(lqtype.eq.1.or.lqtype.eq.2)then
2052         write(6,111)
2053         write(6,110) CHANS1,gproc(3)
2054         write(6,*) 'I        -----------------------        I        I'
2055         write(6,112)
2056         write(6,110) CHANU1,gproc(5)
2057       elseif(lqtype.eq.3)then
2058         write(6,111)
2059         write(6,110) CHANS1,gproc(4)
2060         write(6,*) 'I        -----------------------        I        I'
2061         write(6,112)
2062         write(6,110) CHANU1,gproc(6)
2063       elseif(lqtype.eq.4)then
2064         write(6,111)
2065         write(6,110) CHANS1,gproc(4)
2066         write(6,110) CHANS2,gproc(3)
2067         write(6,*) 'I        -----------------------        I        I'
2068         write(6,112)
2069         write(6,110) CHANU1,gproc(6)
2070         write(6,110) CHANU2,gproc(5)
2071       elseif(lqtype.eq.5)then
2072         write(6,111)
2073         write(6,110) CHANS1,gproc(4)
2074         write(6,*) 'I        -----------------------        I        I'
2075         write(6,112)
2076         write(6,110) CHANU1,gproc(6)
2077       elseif(lqtype.eq.6)then
2078         write(6,111)
2079         write(6,110) CHANS1,gproc(4)
2080         write(6,110) CHANS2,gproc(3)
2081         write(6,*) 'I        -----------------------        I        I'
2082         write(6,112)
2083         write(6,110) CHANU1,gproc(6)
2084         write(6,110) CHANU2,gproc(5)
2085       elseif(lqtype.eq.7)then
2086         write(6,111)
2087         write(6,110) CHANS1,gproc(3)
2088         write(6,*) 'I        -----------------------        I        I'
2089         write(6,112)
2090         write(6,110) CHANU1,gproc(5)
2091       elseif(lqtype.eq.8.or.lqtype.eq.9)then
2092         write(6,111)
2093         write(6,110) CHANS1,gproc(2)
2094         write(6,*) 'I        -----------------------        I        I'
2095         write(6,112)
2096         write(6,110) CHANU1,gproc(8)
2097       elseif(lqtype.eq.10)then
2098         write(6,111)
2099         write(6,110) CHANS1,gproc(1)
2100         write(6,*) 'I        -----------------------        I        I'
2101         write(6,112)
2102         write(6,110) CHANU1,gproc(7)
2103       elseif(lqtype.eq.11)then
2104         write(6,111)
2105         write(6,110) CHANS1,gproc(2)
2106         write(6,110) CHANS2,gproc(1)
2107         write(6,*) 'I        -----------------------        I        I'
2108         write(6,112)
2109         write(6,110) CHANU1,gproc(8)
2110         write(6,110) CHANU2,gproc(7)
2111       elseif(lqtype.eq.12)then
2112         write(6,111)
2113         write(6,110) CHANS1,gproc(1)
2114         write(6,*) 'I        -----------------------        I        I'
2115         write(6,112)
2116         write(6,110) CHANU1,gproc(7)
2117       elseif(lqtype.eq.13)then
2118         write(6,111)
2119         write(6,110) CHANS1,gproc(2)
2120         write(6,110) CHANS2,gproc(1)
2121         write(6,*) 'I        -----------------------        I        I'
2122         write(6,112)
2123         write(6,110) CHANU1,gproc(8)
2124         write(6,110) CHANU2,gproc(7)
2125       elseif(lqtype.eq.14)then
2126         write(6,111)
2127         write(6,110) CHANS1,gproc(2)
2128         write(6,*) 'I        -----------------------        I        I'
2129         write(6,112)
2130         write(6,110) CHANU1,gproc(8)
2131       endif                                                                 
2132       write(6,*) 'I                                       I        I'        
2133       write(6,*) '=================================================='
2134       write(6,*)
2135       write(6,*) '           ---------------------------            '
2136       write(6,*)
2137       write(6,*)
2138       write(6,*) '>>>>>>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<'
2139 101   format('   Number of generated events.............',I8)
2140 C 102   format('   Total cross section (mb)...............',I2)
2141 C 103   format('   s-channel cross section (mb)....... ',F7.2)
2142 104   format('   Number of s-channel generated events...',I8)
2143 C 105   format('   u-channel cross section (mb)....... ',F7.2)
2144 106   format('   Number of u-channel generated events...',I8)
2145 110   format(1X,'I ',A38,'I',1X,I6,' I')
2146 111   format(1X,'I ','s-channel:',28X,'I',7X,' I')
2147 112   format(1X,'I ','u-channel:',28X,'I',7X,' I')
2148       return
2149       end      
2150 *
2151       Subroutine LQGCHAN(CHANS1,CHANS2,CHANU1,CHANU2)
2152 C-----------------------------------------------------
2153 C...Set up the character variables containing 
2154 C...         the generated processes
2155 C-----------------------------------------------------
2156 *
2157       Implicit None
2158 *
2159       Character*38 CHANS1,CHANS2,CHANU1,CHANU2
2160 *
2161 C# LQGpar1.inc #
2162       integer echar
2163       double precision ebeam,pbeam
2164       common /LQGbeam/ebeam,pbeam,echar
2165 
2166 C# LQGproc.inc #
2167       double precision Mlq,G1,G2
2168       Integer LQtype,l_o,q_i,q_j
2169       common /LQGproc/ Mlq,G1,G2,LQtype,l_o,q_i,q_j       
2170 *
2171       Character*7 LQCHA(14)
2172       DATA LQCHA /'S_0L','S_0R','~S_0R','S_1L',
2173      >            'V_1/2L','V_1/2R','~V_1/2L',
2174      >            'V_0L','V_0R','~V_0R','V_1L',
2175      >            'S_1/2L','S_1/2R','~S_1/2L'/                        
2176       character*5 uch(3), dch(3)
2177       character*5 ubch(3), dbch(3)
2178       character*5 lchp(3),lchm(3)
2179       data uch,ubch /' u ',' c ',' t ', ' ubar',' cbar',' tbar'/
2180       data dch,dbch /' d ',' s ',' b ', ' dbar',' sbar',' bbar'/
2181       data lchp,lchm /'e+','mu+','tau+','e-','mu-','tau-'/
2182       character*5 l_in
2183       character*5 q_in
2184       character*5 l_ou
2185       character*5 q_ou
2186 *      
2187       if(echar.gt.0) then
2188        l_in ='e+'
2189        l_ou=lchp(l_o)
2190       else
2191        l_in ='e-'
2192        l_ou=lchm(l_o)
2193       endif
2194       IF(lqtype.eq.1.or.lqtype.eq.2.or.lqtype.eq.7)Then
2195        if(echar.gt.0) then
2196         q_in=ubch(q_i)
2197         q_ou=ubch(q_j)
2198        else
2199         q_in=uch(q_i)
2200         q_ou=uch(q_j)
2201        endif
2202        CHANS1=l_in//q_in//' -> '//LQCHA(lqtype)//' -> '//l_ou//q_ou       
2203        CHANS2=' '
2204        if(echar.gt.0) then
2205         q_in=uch(q_j)
2206         q_ou=uch(q_i)
2207        else
2208         q_in=ubch(q_i)
2209         q_ou=ubch(q_j)
2210        endif
2211        CHANU1=l_in//q_in//' -> '//LQCHA(lqtype)//' -> '//l_ou//q_ou       
2212        CHANU2=' '
2213       ELSEIF(lqtype.eq.3.or.lqtype.eq.5)Then
2214        if(echar.gt.0) then
2215         q_in=dbch(q_i)
2216         q_ou=dbch(q_j)
2217        else
2218         q_in=dch(q_i)
2219         q_ou=dch(q_j)
2220        endif
2221        CHANS1=l_in//q_in//' -> '//LQCHA(lqtype)//' -> '//l_ou//q_ou       
2222        CHANS2=' '
2223        if(echar.gt.0) then
2224         q_in=dch(q_j)
2225         q_ou=dch(q_i)
2226        else
2227         q_in=dbch(q_j)
2228         q_ou=dbch(q_i)        
2229        endif
2230        CHANU1=l_in//q_in//' -> '//LQCHA(lqtype)//' -> '//l_ou//q_ou       
2231        CHANU2=' '
2232       ELSEIF(lqtype.eq.4.or.lqtype.eq.6)Then
2233        if(echar.gt.0) then
2234         q_in=dbch(q_i)
2235         q_ou=dbch(q_j)
2236        else
2237         q_in=dch(q_i)
2238         q_ou=dch(q_j)
2239        endif
2240        CHANS1=l_in//q_in//' -> '//LQCHA(lqtype)//' -> '//l_ou//q_ou       
2241        if(echar.gt.0) then
2242         q_in=ubch(q_i)
2243         q_ou=ubch(q_j)
2244        else
2245         q_in=uch(q_i)
2246         q_ou=uch(q_j)
2247        endif
2248        CHANS2=l_in//q_in//' -> '//LQCHA(lqtype)//' -> '//l_ou//q_ou
2249        if(echar.gt.0) then
2250         q_in=dch(q_j)
2251         q_ou=dch(q_i)
2252        else
2253         q_in=dbch(q_j)
2254         q_ou=dbch(q_i)
2255        endif
2256        CHANU1=l_in//q_in//' -> '//LQCHA(lqtype)//' -> '//l_ou//q_ou       
2257        if(echar.gt.0) then
2258         q_in=uch(q_j)
2259         q_ou=uch(q_i)
2260        else
2261         q_in=ubch(q_j)
2262         q_ou=ubch(q_i)
2263        endif
2264        CHANU2=l_in//q_in//' -> '//LQCHA(lqtype)//' -> '//l_ou//q_ou       
2265       ELSEIF(lqtype.eq.8.or.lqtype.eq.9.or.lqtype.eq.14)Then
2266        if(echar.gt.0) then
2267         q_in=dch(q_i)
2268         q_ou=dch(q_j)
2269        else
2270         q_in=dbch(q_i)
2271         q_ou=dbch(q_j)
2272        endif
2273        CHANS1=l_in//q_in//' -> '//LQCHA(lqtype)//' -> '//l_ou//q_ou       
2274        CHANS2=' '
2275        if(echar.gt.0) then
2276         q_in=dbch(q_j)
2277         q_ou=dbch(q_i)
2278        else
2279         q_in=dch(q_i)
2280         q_ou=dch(q_j)
2281        endif
2282        CHANU1=l_in//q_in//' -> '//LQCHA(lqtype)//' -> '//l_ou//q_ou       
2283        CHANU2=' '
2284       ELSEIF(lqtype.eq.10.or.lqtype.eq.12)Then
2285        if(echar.gt.0) then
2286         q_in=uch(q_i)
2287         q_ou=uch(q_j)
2288        else
2289         q_in=ubch(q_i)
2290         q_ou=ubch(q_j)
2291        endif
2292        CHANS1=l_in//q_in//' -> '//LQCHA(lqtype)//' -> '//l_ou//q_ou       
2293        CHANS2=' '
2294        if(echar.gt.0) then
2295         q_in=ubch(q_j)
2296         q_ou=ubch(q_i)
2297        else
2298         q_in=uch(q_i)
2299         q_ou=uch(q_j)
2300        endif
2301        CHANU1=l_in//q_in//' -> '//LQCHA(lqtype)//' -> '//l_ou//q_ou       
2302        CHANU2=' '       
2303       ELSEIF(lqtype.eq.11.or.lqtype.eq.13)Then
2304        if(echar.gt.0) then
2305         q_in=dch(q_i)
2306         q_ou=dch(q_j)
2307        else
2308         q_in=dbch(q_i)
2309         q_ou=dbch(q_j)
2310        endif
2311        CHANS1=l_in//q_in//' -> '//LQCHA(lqtype)//' -> '//l_ou//q_ou       
2312        if(echar.gt.0) then
2313         q_in=uch(q_i)
2314         q_ou=uch(q_j)
2315        else
2316         q_in=ubch(q_i)
2317         q_ou=ubch(q_j)
2318        endif
2319        CHANS2=l_in//q_in//' -> '//LQCHA(lqtype)//' -> '//l_ou//q_ou       
2320        if(echar.gt.0) then
2321         q_in=dbch(q_j)
2322         q_ou=dbch(q_i)
2323        else
2324         q_in=dch(q_j)
2325         q_ou=dch(q_i)
2326        endif
2327        CHANU1=l_in//q_in//' -> '//LQCHA(lqtype)//' -> '//l_ou//q_ou       
2328        if(echar.gt.0) then
2329         q_in=ubch(q_j)
2330         q_ou=ubch(q_i)
2331        else
2332         q_in=uch(q_j)
2333         q_ou=uch(q_i)
2334        endif
2335        CHANU2=l_in//q_in//' -> '//LQCHA(lqtype)//' -> '//l_ou//q_ou       
2336       ENDIF      
2337       Return
2338       End
2339 *
2340 cc      Subroutine LQGTESTL
2341 C-------------------------------------------
2342 C...   Test routine to generate the process
2343 C...   e q_1 -> S_1/2^L -> mu q_1
2344 C...   Mass of the LQ = 250 GeV
2345 C...   beams at HERA energies
2346 C-------------------------------------------
2347 ccc      Implicit none 
2348 C...Pythia parameters
2349 C...Parameters. 
2350 cc      double precision PARP,PARI
2351 cc      integer MSTP,MSTI
2352 cc      COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
2353 
2354 C...LQGENEP run setup parameters
2355 cc       double precision BEAMPAR,LQGPAR3
2356 cc       integer LQGPAR1,LQGPAR2
2357 cc       COMMON/LQGDATA/BEAMPAR(3),LQGPAR1(10),LQGPAR2(10),LQGPAR3(20)
2358 
2359 
2360 C...LQGENEP event informations
2361 cc       double precision LQGKPAR,LQGDDX
2362 cc       integer LQGPID
2363 cc       COMMON/LQGEVT/LQGKPAR(3),LQGDDX(3),LQGPID(3)
2364 *
2365 
2366 cc       integer NeV,i
2367 *
2368 cc       beampar(1)=1.
2369 cc       beampar(2)=27.5
2370 cc       beampar(3)=820.
2371 cc       lqgpar1(1)=5000
2372 cc       lqgpar1(2)=10
2373 cc       lqgpar1(3)=1
2374 cc       lqgpar1(4)=0
2375 cc       lqgpar1(5)=0
2376 cc       lqgpar2(1)=12
2377 cc       lqgpar2(2)=1
2378 cc       lqgpar2(3)=1
2379 cc       lqgpar2(4)=2
2380 cc       lqgpar3(1)=250.
2381 cc       lqgpar3(2)=0.3
2382 cc       lqgpar3(3)=0.3
2383 cc       lqgpar3(4)=0.
2384 cc       lqgpar3(5)=1.
2385 cc       lqgpar3(6)=0.
2386 cc       lqgpar3(7)=1.cc       
2387 cc       lqgpar3(8)=500.
2388 cc       lqgpar3(9)=1.cc       
2389 cc       lqgpar3(10)=4.cc       
2390 cc       lqgpar3(11)=32.cc       
2391 * Max cross section
2392 cc       lqgpar3(12)=3.d-6
2393 *
2394 * switch off initial state QCD and QED radiation
2395 ccc       MSTP(61)=0
2396 * switch off final state QCD and QED radiation
2397 ccc       MSTP(71)=0
2398 * switch off multiple interaction
2399 ccc       MSTP(81)=0
2400 * switch off fragmentation and decay
2401 ccc       MSTP(111)=0
2402 
2403 * LQGENEP Initialization
2404 cc       call LQGENEP(0)
2405 cc       Nev=lqgpar1(1)
2406 
2407 * LQGENEP generation loop
2408 cc       do i=1,Nev
2409 cc        call LQGENEP(1)
2410 cc       enddo
2411 
2412 * LQGENEP termination
2413 cc       call LQGENEP(2)
2414 
2415 cc       return
2416 cc       end
2417 *
2418       Subroutine LQGTESTH
2419 ***************************
2420 C-------------------------------------------
2421 C...   Test routine to generate the process
2422 C...   e+ q_2 -> ~S_0^R -> mu+ q_1
2423 C...   Mass of the LQ = 600 GeV
2424 C...   beams at HERA energies
2425 C-------------------------------------------
2426       implicit none 
2427 C...Pythia parameters
2428 C...Parameters. 
2429       double precision PARP,PARI
2430       integer MSTP,MSTI
2431       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
2432 
2433 C...LQGENEP run setup parameters
2434        double precision BEAMPAR,LQGPAR3
2435        integer LQGPAR1,LQGPAR2
2436        CHARACTER*256 out_file
2437        COMMON/LQGDATA/BEAMPAR(3),LQGPAR1(10),LQGPAR2(10),LQGPAR3(20),
2438      > out_file
2439 
2440 
2441 
2442 C...LQGENEP event informations
2443        double precision LQGKPAR,LQGDDX
2444        integer LQGPID
2445        COMMON/LQGEVT/LQGKPAR(3),LQGDDX(3),LQGPID(3)
2446 *
2447        integer NeV,i
2448 *
2449        beampar(1)=1.
2450        beampar(2)=27.5
2451        beampar(3)=820.
2452        lqgpar1(1)=5000
2453        lqgpar1(2)=10
2454        lqgpar1(3)=1
2455        lqgpar1(4)=0
2456        lqgpar1(5)=0
2457        lqgpar2(1)=3
2458        lqgpar2(2)=2
2459        lqgpar2(3)=1
2460        lqgpar2(4)=3
2461        lqgpar3(1)=600.
2462        lqgpar3(2)=0.3
2463        lqgpar3(3)=0.3
2464        lqgpar3(4)=0.
2465        lqgpar3(5)=1.
2466        lqgpar3(6)=0.
2467        lqgpar3(7)=1.       
2468        lqgpar3(8)=500.
2469        lqgpar3(9)=1.       
2470        lqgpar3(10)=4.       
2471        lqgpar3(11)=32.       
2472 * Max cross section
2473        lqgpar3(12)=2.d-11
2474 *
2475 * switch off initial state QCD and QED radiation
2476 c       MSTP(61)=0
2477 * switch off final state QCD and QED radiation
2478 c       MSTP(71)=0
2479 * switch off multiple interaction
2480 c       MSTP(81)=0
2481 * switch off fragmentation and decay
2482 c       MSTP(111)=0
2483 
2484 * LQGENEP Initialization
2485 cc       call LQGENEP(0)
2486 cc       Nev=lqgpar1(1)
2487 
2488 * LQGENEP generation loop
2489 cc       do i=1,Nev
2490 cc        call LQGENEP(1)
2491 cc       enddo       
2492 
2493 * LQGENEP termination
2494 cc       call LQGENEP(2)
2495 
2496        return
2497        end
2498 *
2499       Subroutine LQGXNWA(XNWA,IERR)
2500 C-------------------------------------------------------      
2501 C...Evaluates Narrow Width Approximation Cross Section  
2502 C...                                                      
2503 C...  output argument: XNWA = NWA cross section           
2504 C...                   IERR = 0 -> OK , 1 -> LQ higher 
2505 C...                                    than center of 
2506 C...                                    mass energy 
2507 C-------------------------------------------------------      
2508       Implicit None
2509 *
2510       double precision XNWA
2511       integer IERR
2512       double precision xx,cutfac0,cutfac1,ycut,factor
2513       double precision upv,dnv,usea,dsea,str,chm,bot,top,gl
2514       double precision xsu,xsd,xsub,xsdb
2515       double precision Br_rat
2516       
2517 C# LQGpar1.inc #
2518       integer echar
2519       double precision ebeam,pbeam
2520       common /LQGbeam/ebeam,pbeam,echar
2521 
2522 C# LQGKinC.inc #
2523       double precision xmax,xmin,ymax,ymin,zmax,zmin,Q2min
2524       common /LQGKinC/ xmax,xmin,ymax,ymin,zmax,zmin,Q2min
2525 
2526 C# LQGproc.inc #
2527       double precision Mlq,G1,G2
2528       Integer LQtype,l_o,q_i,q_j
2529       common /LQGproc/ Mlq,G1,G2,LQtype,l_o,q_i,q_j 
2530 
2531 C# LQGKinV.inc # 
2532       double precision S,Srad,x,y,z,Q2
2533       common /LQGKinV/ S,Srad,x,y,z,Q2
2534 *
2535       double precision pi
2536       parameter (pi=3.141592653589793d0)
2537 *
2538       IERR=0
2539       XNWA=0.d0
2540       xx=MLQ**2/S
2541       if (xx.gt.1d0.or.xx.lt.1d-6) then
2542        write(6,*)'#################################################'
2543        write(6,*)'LQGXNWA Error:'
2544        write(6,*)'LQ mass not compatible with center of mass energy'
2545        write(6,*)'#################################################'
2546        IERR=1
2547        return
2548       endif
2549       if (Q2min.le.1.) then 
2550        cutfac0=1.d0
2551        cutfac1=2.d0
2552       else
2553        ycut=Q2min/xx/s
2554        if (ycut.ge.1d0) then
2555         XNWA=0.d0
2556        else
2557         cutfac0=1.d0-ycut
2558         cutfac1=2.*(cutfac0)**3
2559        endif
2560       endif
2561       factor=(pi*G1**2)/(4.d0*s)*(0.3894)
2562       call structm(xx,mlq,upv,dnv,usea,dsea,str,chm,bot,top,gl)
2563       if (echar.ge.0) then
2564        xsu=factor*(upv+usea)/xx
2565        xsd=factor*(dnv+dsea)/xx
2566        xsub=factor*usea/xx
2567        xsdb=factor*dsea/xx
2568       else
2569        xsu=factor*usea/xx
2570        xsd=factor*dsea/xx
2571        xsub=factor*(upv+usea)/xx
2572        xsdb=factor*(dnv+dsea)/xx
2573       endif
2574       IF(l_o.eq.1)then
2575 * only electron and eventually neutrino final states considered      
2576        IF(lqtype.eq.1)XNWA = cutfac0*xsub*0.5
2577        IF(lqtype.eq.2)XNWA = cutfac0*xsub
2578        IF(lqtype.eq.3)XNWA = cutfac0*xsdb
2579        IF(lqtype.eq.4)XNWA = cutfac0*(xsub*0.5+2.*xsdb)
2580        IF(lqtype.eq.12)XNWA = cutfac0*xsu
2581        IF(lqtype.eq.13)XNWA = cutfac0*(xsu+xsd)
2582        IF(lqtype.eq.14)XNWA = cutfac0*xsd
2583        IF(lqtype.eq.5)XNWA = cutfac1*(xsdb)
2584        IF(lqtype.eq.6)XNWA = cutfac1*(xsub+xsdb)
2585        IF(lqtype.eq.7)XNWA = cutfac1*(xsub)*0.5
2586        IF(lqtype.eq.8)XNWA = cutfac1*(xsd)*0.5
2587        IF(lqtype.eq.9)XNWA = cutfac1*(xsd)
2588        IF(lqtype.eq.10)XNWA = cutfac1*(xsu)
2589        IF(lqtype.eq.11)XNWA = cutfac1*(2.*xsu+xsd*0.5)
2590       ELSE
2591 * electron and muon (or tau) possible final states
2592        Br_rat=G2**2/(G1**2+G2**2)      
2593        IF(lqtype.eq.1)XNWA = Br_rat*cutfac0*xsub*0.5
2594        IF(lqtype.eq.2)XNWA = Br_rat*cutfac0*xsub
2595        IF(lqtype.eq.3)XNWA = Br_rat*cutfac0*xsdb
2596        IF(lqtype.eq.4)XNWA = Br_rat*cutfac0*(xsub*0.5+2.*xsdb)
2597        IF(lqtype.eq.12)XNWA = Br_rat*cutfac0*xsu
2598        IF(lqtype.eq.13)XNWA = Br_rat*cutfac0*(xsu+xsd)
2599        IF(lqtype.eq.14)XNWA = Br_rat*cutfac0*xsd
2600        IF(lqtype.eq.5)XNWA = Br_rat*cutfac1*(xsdb)
2601        IF(lqtype.eq.6)XNWA = Br_rat*cutfac1*(xsub+xsdb)
2602        IF(lqtype.eq.7)XNWA = Br_rat*cutfac1*(xsub)*0.5
2603        IF(lqtype.eq.8)XNWA = Br_rat*cutfac1*(xsd)*0.5
2604        IF(lqtype.eq.9)XNWA = Br_rat*cutfac1*(xsd)
2605        IF(lqtype.eq.10)XNWA = Br_rat*cutfac1*(xsu)
2606        IF(lqtype.eq.11)XNWA = Br_rat*cutfac1*(2.*xsu+xsd*0.5)
2607       ENDIF
2608       return
2609       end
2610 *
2611       subroutine LQGXINT(XINT,XINTE,IERR)
2612 C-------------------------------------------------------      
2613 C... Cross Section evaluation by Double Differential      
2614 C... Cross Section integration                            
2615 C...                                                      
2616 C... Output argument: XINT =  Integrated Cross Section    
2617 C...                  XINTE = Error on Cross Section      
2618 C...                  IERR =  0 -> OK , >0 -> problems in  
2619 C...                            cross section evaluation    
2620 C-------------------------------------------------------      
2621       Implicit None
2622 *
2623       Double precision XINT,XINTE
2624       Integer IERR
2625 *
2626       double precision relerr,result
2627       external DSDXDY
2628       integer minpts,maxpts,iwk,nfnevl,ifail
2629       parameter (minpts=10000)
2630       parameter (maxpts=5000000)
2631       parameter (iwk=1000000)
2632       double precision wk(iwk),a(2),b(2),eps
2633       parameter (eps=1d-3) 
2634 
2635 C# LQGKinC.inc #
2636       double precision xmax,xmin,ymax,ymin,zmax,zmin,Q2min
2637       common /LQGKinC/ xmax,xmin,ymax,ymin,zmax,zmin,Q2min
2638 
2639       IERR=0
2640       XINT=0.d0
2641       XINTE=0.d0
2642 *
2643       a(1)=XMIN
2644       b(1)=XMAX
2645       a(2)=YMIN
2646       b(2)=YMAX
2647 *      
2648       call dadmul(DSDXDY,2,a,b,minpts,maxpts,eps,
2649      +               wk,iwk,result,relerr,nfnevl,ifail)
2650 *
2651       if(ifail.ne.0)then
2652        Write(6,*)'LQGXINT: Error in Cross Section integration'
2653        IERR=ifail
2654        return
2655       endif     
2656 *
2657       XINT=result
2658       XINTE=abs(result)*relerr
2659 *
2660       return
2661       end
2662 *
2663       Double precision function dsdxdy(n,xx)
2664 *
2665       IMPLICIT None
2666       
2667       integer n
2668       double precision xx(*)
2669 
2670 C# LQGKinC.inc #
2671       double precision xmax,xmin,ymax,ymin,zmax,zmin,Q2min
2672       common /LQGKinC/ xmax,xmin,ymax,ymin,zmax,zmin,Q2min
2673 
2674 C# LQGproc.inc #
2675       double precision Mlq,G1,G2
2676       Integer LQtype,l_o,q_i,q_j
2677       common /LQGproc/ Mlq,G1,G2,LQtype,l_o,q_i,q_j 
2678 
2679 C# LQGKinV.inc # 
2680       double precision S,Srad,x,y,z,Q2
2681       common /LQGKinV/ S,Srad,x,y,z,Q2
2682 
2683 C# LQGout.inc #
2684       double precision DXSEC(3),pvalence
2685       integer q_o,q_s,genproc,genprtyp,sch,uch,gproc(8)
2686       common /LQGout/ DXSEC,pvalence,q_o,q_s,genproc,genprtyp,
2687      >sch,uch,gproc
2688 
2689 * conversion from pb to mb
2690       double precision CONV
2691       DATA CONV/1.D-9/
2692 *
2693       DSDXDY=0.
2694       X=xx(1)
2695       Y=xx(2)
2696       Q2=x*y*s        
2697       dsdxdy=0.d0
2698       if(Q2.gt.Q2min) then             
2699        call LQGDDXS
2700        DSDXDY=(DXSEC(2)+DXSEC(3))*conv
2701       endif
2702       return
2703       end
2704 *
2705       Subroutine LQGXHMA(XSEC,XSECE,IERR)
2706 C---------------------------------------------------------      
2707 C...    Evaluates High Mass Approximation Cross Section     
2708 C...                                                        
2709 C...  Output arguments: XSEC  = HMA cross section (mb)      
2710 C...                    XSECE = error on cross section (mb) 
2711 C...                    IERR = 0 -> OK , >0 -> problems in   
2712 C...                             cross section evaluation      
2713 C...                                                        
2714 C---------------------------------------------------------      
2715 *
2716       implicit none
2717 *     
2718       double precision XSEC,XSECE
2719       integer IERR
2720       external qdens
2721 
2722 * declaration for DADMUL variables
2723       integer minpts,maxpts,iwk
2724       parameter (minpts=10000)
2725       parameter (maxpts=5000000)
2726       parameter (iwk=1000000)
2727       double precision wk(iwk),a(2),b(2)
2728       double precision eps
2729       parameter (eps=1d-3)  
2730       double precision result,relerr
2731       integer nfnevl,ifail
2732       double precision pi
2733       parameter (pi=3.141592653589793d0)
2734 *
2735 C# LQGKinC.inc #
2736       double precision xmax,xmin,ymax,ymin,zmax,zmin,Q2min
2737       common /LQGKinC/ xmax,xmin,ymax,ymin,zmax,zmin,Q2min
2738 
2739 C# LQGproc.inc #
2740       double precision Mlq,G1,G2
2741       Integer LQtype,l_o,q_i,q_j
2742       common /LQGproc/ Mlq,G1,G2,LQtype,l_o,q_i,q_j 
2743 
2744 C# LQGKinV.inc # 
2745       double precision S,Srad,x,y,z,Q2
2746       common /LQGKinV/ S,Srad,x,y,z,Q2
2747 
2748 *
2749       IERR=0
2750       XSEC=0.d0
2751       XSECE=0.d0
2752       a(1)=xmin
2753       b(1)=xmax
2754       a(2)=ymin
2755       b(2)=ymax
2756 *      
2757       call dadmul(qdens,2,a,b,minpts,maxpts,eps,
2758      +            wk,iwk,result,relerr,nfnevl,ifail)         
2759 *
2760       if(ifail.ne.0)then
2761        Write(6,*)'LQGXHMA: Error in Cross Section integration'
2762        IERR=ifail
2763        return
2764       endif     
2765       XSEC=S/32./pi*G1*G1*G2*G2/MLQ/MLQ/MLQ/MLQ*result*0.38938
2766       XSECE=relerr*XSEC
2767       return
2768       end
2769 *
2770       double precision function qdens(n,xx)
2771 *
2772       implicit none
2773 *
2774       integer n
2775       double precision xx(*)
2776       double precision UPVs,DNVs,USEAs,DSEAs,STRs,
2777      +CHMs,BOTs,TOPs,GLs
2778       double precision UPVu,DNVu,USEAu,DSEAu,STRu,
2779      +CHMu,BOTu,TOPu,GLu
2780 *
2781       double precision sh,u
2782       double precision UPQs(3),DNQs(3),UPQBs(3),DNQBs(3)
2783       double precision UPQu(3),DNQu(3),UPQBu(3),DNQBu(3)
2784 *
2785 C# LQGKinC.inc #
2786       double precision xmax,xmin,ymax,ymin,zmax,zmin,Q2min
2787       common /LQGKinC/ xmax,xmin,ymax,ymin,zmax,zmin,Q2min
2788 
2789 C# LQGproc.inc #
2790       double precision Mlq,G1,G2
2791       Integer LQtype,l_o,q_i,q_j
2792       common /LQGproc/ Mlq,G1,G2,LQtype,l_o,q_i,q_j 
2793 
2794 C# LQGKinV.inc # 
2795       double precision S,Srad,x,y,z,Q2
2796       common /LQGKinV/ S,Srad,x,y,z,Q2
2797 
2798 C# LQGpar1.inc #
2799       integer echar
2800       double precision ebeam,pbeam
2801       common /LQGbeam/ebeam,pbeam,echar
2802 *
2803       X=xx(1)
2804       Y=xx(2)
2805 *      
2806       qdens=0.d0
2807       Q2=S*x*y
2808       sh=sqrt(S*X)
2809       u=sqrt(S*X*(1-Y))
2810       if(Q2.gt.Q2min)Then
2811        CALL STRUCTM(X,sh,UPVs
2812      &            ,DNVs,USEAs,DSEAs,STRs,CHMs,BOTs,TOPs,GLs) 
2813       else
2814        UPVs=0.
2815        DNVs=0.
2816        USEAs=0.
2817        DSEAs=0.
2818        STRs=0.
2819        CHMs=0.
2820        BOTs=0.
2821        TOPs=0.
2822        GLs=0.        
2823       endif
2824       if(Q2.gt.Q2min)Then
2825        CALL STRUCTM(X,u,UPVu
2826      &            ,DNVu,USEAu,DSEAu,STRu,CHMu,BOTu,TOPu,GLu)
2827       else
2828        UPVu=0.
2829        DNVu=0.
2830        USEAu=0.
2831        DSEAu=0.
2832        STRu=0.
2833        CHMu=0.
2834        BOTu=0.
2835        TOPu=0.
2836        GLu=0.        
2837       endif 
2838 *
2839       if(echar.eq.1)then
2840 * case e+, mu+, tau+ 
2841        UPQu(1)=UPVu+USEAu
2842        UPQBu(1)=USEAu
2843        UPQu(2)=CHMu
2844        UPQBu(2)=CHMu
2845        UPQu(3)=TOPu
2846        UPQBu(3)=TOPu
2847        DNQu(1)=DNVu+DSEAu
2848        DNQBu(1)=DSEAu
2849        DNQu(2)=STRu
2850        DNQBu(2)=STRu
2851        DNQu(3)=BOTu
2852        DNQBu(3)=BOTu
2853       elseif(echar.eq.-1)then
2854 * case e-, mu-, tau- 
2855        UPQu(1)=USEAu
2856        UPQBu(1)=UPVu+USEAu
2857        UPQu(2)=CHMu
2858        UPQBu(2)=CHMu
2859        UPQu(3)=TOPu
2860        UPQBu(3)=TOPu
2861        DNQu(1)=DSEAu
2862        DNQBu(1)=DNVu+DSEAu
2863        DNQu(2)=STRu
2864        DNQBu(2)=STRu
2865        DNQu(3)=BOTu
2866        DNQBu(3)=BOTu
2867       endif          
2868 * s channel densities
2869       if(echar.eq.1)then
2870 * case e+, mu+, tau+ 
2871        UPQs(1)=UPVs+USEAs
2872        UPQBs(1)=USEAs
2873        UPQs(2)=CHMs
2874        UPQBs(2)=CHMs
2875        UPQs(3)=TOPs
2876        UPQBs(3)=TOPs
2877        DNQs(1)=DNVs+DSEAs
2878        DNQBs(1)=DSEAs
2879        DNQs(2)=STRs
2880        DNQBs(2)=STRs
2881        DNQs(3)=BOTs
2882        DNQBs(3)=BOTs
2883       elseif(echar.eq.-1)then
2884 * case e-, mu-, tau- 
2885        UPQs(1)=USEAs
2886        UPQBs(1)=UPVs+USEAs
2887        UPQs(2)=CHMs
2888        UPQBs(2)=CHMs
2889        UPQs(3)=TOPs
2890        UPQBs(3)=TOPs
2891        DNQs(1)=DSEAs
2892        DNQBs(1)=DNVs+DSEAs
2893        DNQs(2)=STRs
2894        DNQBs(2)=STRs
2895        DNQs(3)=BOTs
2896        DNQBs(3)=BOTs
2897       endif
2898 *
2899       if(LQTYPE.eq.1)Then
2900        qdens=UPQBs(q_i)/2.+UPQu(q_j)*(1-y)*(1-y)/2.
2901       elseif(LQTYPE.eq.2)Then
2902        qdens=UPQBs(q_i)/2.+UPQu(q_j)*(1-y)*(1-y)/2.
2903       elseif(LQTYPE.eq.3)Then
2904        qdens=DNQBs(q_i)/2.+DNQu(q_j)*(1-y)*(1-y)/2.
2905       elseif(LQTYPE.eq.4)Then
2906        if(q_j.eq.3)UPQBs(q_i)=0.
2907        if(q_i.eq.3)UPQu(q_j)=0.
2908        qdens=(UPQBs(q_i)+4.*DNQBs(q_i))/2.
2909      >       +(UPQu(q_j)+4.*DNQu(q_j))*(1.-Y)*(1.-Y)/2.
2910       elseif(LQTYPE.eq.5)Then
2911        qdens=DNQBs(q_i)*(1-y)*(1-y)*2.+DNQu(q_j)*2.
2912       elseif(LQTYPE.eq.6)Then
2913        if(q_j.eq.3)UPQBs(q_i)=0.
2914        if(q_i.eq.3)UPQu(q_j)=0.
2915        qdens=(UPQBs(q_i)+DNQBs(q_i))*(1-y)*(1-y)*2.
2916      >       +(UPQu(q_j)+DNQu(q_j))*2.
2917       elseif(LQTYPE.eq.7)Then
2918        qdens=UPQBs(q_i)*(1-y)*(1-y)*2.+UPQu(q_j)*2.
2919       elseif(LQTYPE.eq.8)Then
2920        qdens=DNQs(q_i)*2.*(1.-Y)*(1.-Y)+DNQBu(q_j)*2.
2921       elseif(LQTYPE.eq.9)Then
2922        qdens=DNQs(q_i)*2.*(1.-Y)*(1.-Y)+DNQBu(q_j)*2.
2923       elseif(LQTYPE.eq.10)Then  
2924        qdens=UPQs(q_i)*2.*(1.-Y)*(1.-Y)+UPQBu(q_j)*2.
2925       elseif(LQTYPE.eq.11)Then  
2926        if(q_j.eq.3)UPQs(q_i)=0.
2927        if(q_i.eq.3)UPQBu(q_j)=0.
2928        qdens=(4.*UPQs(q_i)+DNQs(q_i))*2.*(1.-Y)*(1.-Y)
2929      >      +(4.*UPQBu(q_j)+DNQBu(q_j))*2.
2930       elseif(LQTYPE.EQ.12)Then
2931        qdens=UPQs(q_i)/2.+UPQBu(q_j)*(1.-Y)*(1.-Y)/2.
2932       elseif(LQTYPE.EQ.13)Then
2933        if(q_j.eq.3)UPQs(q_i)=0.
2934        if(q_i.eq.3)UPQBu(q_j)=0.
2935        qdens=(UPQs(q_i)+DNQs(q_i))/2.
2936      >       +(UPQBu(q_j)+DNQBu(q_j))*(1.-Y)*(1.-Y)/2.
2937       elseif(LQTYPE.EQ.14)Then
2938        qdens=DNQs(q_i)/2.+DNQBu(q_j)*(1.-Y)*(1.-Y)/2.
2939       endif
2940       return
2941       end
2942 *