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