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