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