source/LQGENEP/lqgenep_EIC.f

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