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