Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-17 09:57:06 

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

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!