67 common/
pydat1/mstu(200),paru(200),mstj(200),parj(200)
69 DOUBLE PRECISION paru,parj
70 common/
pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
71 INTEGER mdcy,mdme,kfdp
73 common/
pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
74 INTEGER msel,mselpd,msub,kfin
79 common/
pydatr/mrpy(6),rrpy(100)
83 common/hepmcid/hpmcfid,logfid
84 integer hpmcfid,logfid
86 common/npdf/
mass,nset,eps09,initstr
95 common/evrecord/nsim,
npart,
offset,hadrotype,sqrts,collider,hadro,
98 double precision sqrts
100 character*2 isochannel
101 logical hadro,shorthepmc
103 common/disc/ndisc,nstrange,ngood,errcount,wdisc,discard
105 INTEGER ndisc,nstrange,ngood,errcount
106 double precision wdisc
108 common/
weight/evweight,sumofweights
109 double precision evweight,sumofweights
111 common/check/nscat,nscateff,nsplit
112 DOUBLE PRECISION nscat,nscateff,nsplit
114 common/extrapolations/ntotspliti,noverspliti,ntotpdf,noverpdf,
115 &ntotxsec,noverxsec,ntotsuda,noversuda
116 integer ntotspliti,noverspliti,ntotpdf,noverpdf,
117 &ntotxsec,noverxsec,ntotsuda,noversuda
120 integer nsimpp,nsimpn,nsimnp,nsimnn,nsimsum,nsimchn
121 double precision sumofweightstot,wdisctot,scalefac
131 if (collider.eq.
'EEJJ')
then
136 &
'####################################################'
138 write(logfid,*)
'generating ',nsim,
' events in ',
b1,
' + ',b2,
142 &
'####################################################'
152 sumofweightstot = sumofweightstot+sumofweights
153 wdisctot = wdisctot + wdisc
155 write(logfid,*)
'cross section in e+ + e- channel:',
pari(1),
'mb'
156 write(logfid,*)
'sum of event weights in e+ + e- channel:',
162 if (isochannel.eq.
'PP')
then
167 elseif (isochannel.eq.
'PN')
then
172 elseif (isochannel.eq.
'NP')
then
177 elseif (isochannel.eq.
'NN')
then
187 nsimsum = nsimpp + nsimpn + nsimnp + nsimnn
188 scalefac = nsim*1.d0/(nsimsum*1.d0)
189 nsimpp = int(nsimpp*scalefac)
190 nsimpn = int(nsimpn*scalefac)
191 nsimnp = int(nsimnp*scalefac)
192 nsimnn = int(nsimnn*scalefac)
193 nsimsum = nsimpp + nsimpn + nsimnp + nsimnn
201 elseif (kk.eq.2)
then
205 elseif (kk.eq.3)
then
216 &
'####################################################'
218 write(logfid,*)
'generating ',nsimchn,
' events in ',
219 &
b1,
' + ',b2,
' channel'
222 &
'####################################################'
232 sumofweightstot = sumofweightstot+sumofweights
233 wdisctot = wdisctot + wdisc
235 write(logfid,*)
'cross section in ',
b1,
' + ',b2,
' channel:',
237 write(logfid,*)
'sum of event weights in ',
b1,
' + ',b2,
238 &
' channel:',sumofweights-wdisc
244 WRITE(hpmcfid,
'(A)')
'HepMC::IO_GenEvent-END_EVENT_LISTING'
246 CLOSE(hpmcfid,status=
'keep')
249 write(logfid,*)
'mean number of scatterings:',
250 & nscat/(sumofweightstot-wdisctot)
251 write(logfid,*)
'mean number of effective scatterings:',
252 & nscateff/(sumofweightstot-wdisctot)
253 write(logfid,*)
'mean number of splittings:',
254 & nsplit/(sumofweightstot-wdisctot)
256 write(logfid,*)
'number of extrapolations in splitting integral: ',
257 & noverspliti,
' (',(noverspliti*1.d0)/(ntotspliti*1.d0),
'%)'
259 &
'number of extrapolations in splitting partonic PDFs: ',
260 & noverpdf,
' (',(noverpdf*1.d0)/(ntotpdf*1.d0),
'%)'
262 &
'number of extrapolations in splitting cross sections: ',
263 & noverxsec,
' (',(noverxsec*1.d0)/(ntotxsec*1.d0),
'%)'
265 &
'number of extrapolations in Sudakov form factor: ',
266 & noversuda,
' (',(noversuda*1.d0)/(ntotsuda*1.d0),
'%)'
268 write(logfid,*)
'number of good events: ',ngood
269 write(logfid,*)
'total number of discarded events: ',ndisc
270 write(logfid,*)
'number of events for which conversion '//
271 &
'to hepmc failed: ',nstrange
274 close(logfid,status=
'keep')
298 common/
pydat1/mstu(200),paru(200),mstj(200),parj(200)
300 DOUBLE PRECISION paru,parj
301 common/
pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
302 INTEGER mdcy,mdme,kfdp
303 DOUBLE PRECISION brat
304 common/
pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
305 INTEGER msel,mselpd,msub,kfin
306 DOUBLE PRECISION ckin
310 common/
pydatr/mrpy(6),rrpy(100)
312 DOUBLE PRECISION rrpy
314 common/npdf/
mass,nset,eps09,initstr
316 DOUBLE PRECISION mass
326 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
329 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
330 LOGICAL angord,scatrecoil,allhad,
compress
332 common/splitint/splitiggv(1000,1000),splitiqqv(1000,1000),
333 &splitiqgv(1000,1000),qval(1000),zmval(1000),qmax,zmmin,npoint
335 DOUBLE PRECISION splitiggv,splitiqqv,splitiqgv,
336 &qval,zmval,qmax,zmmin
338 common/pdfs/qinqx(2,1000),ginqx(2,1000),qingx(2,1000),
340 DOUBLE PRECISION qinqx,ginqx,qingx,gingx
342 common/xsecs/intq1(1001,101),intq2(1001,101),
343 &intg1(1001,101),intg2(1001,101)
344 DOUBLE PRECISION intq1,intq2,intg1,intg2
346 common/insuda/sudaqq(1000,2),sudaqg(1000,2),sudagg(1000,2)
348 DOUBLE PRECISION sudaqq,sudaqg,sudagg,sudagc
350 common/expint/
eix(3,1000),valmax,nval
352 DOUBLE PRECISION eix,valmax
354 common/disc/ndisc,nstrange,ngood,errcount,wdisc,discard
356 INTEGER ndisc,nstrange,ngood,errcount
357 double precision wdisc
359 common/ftimefac/ftfac
360 DOUBLE PRECISION ftfac
362 common/alphasfac/ptfac
363 DOUBLE PRECISION ptfac
365 common/check/nscat,nscateff,nsplit
366 DOUBLE PRECISION nscat,nscateff,nsplit
368 common/extrapolations/ntotspliti,noverspliti,ntotpdf,noverpdf,
369 &ntotxsec,noverxsec,ntotsuda,noversuda
370 integer ntotspliti,noverspliti,ntotpdf,noverpdf,
371 &ntotxsec,noverxsec,ntotsuda,noversuda
373 common/
weight/evweight,sumofweights
374 double precision evweight,sumofweights
376 common/wexpo/weightex
377 DOUBLE PRECISION weightex
379 common/hepmcid/hpmcfid,logfid
380 integer hpmcfid,logfid
388 common/evrecord/nsim,
npart,
offset,hadrotype,sqrts,collider,hadro,
391 double precision sqrts
393 character*2 isochannel
394 logical hadro,shorthepmc
396 common/storescatcen/nscatcen,maxnscatcen,scatflav(10000),
397 &scatcen(10000,5),writescatcen,writedummies
398 integer nscatcen,maxnscatcen,scatflav
399 double precision scatcen
400 logical writescatcen,writedummies
402 common/pythiaparams/ptmin,ptmax,weighted
403 double precision ptmin,ptmax
407 INTEGER njob,ios,
pos,
i,
j,jj,intmass
411 CHARACTER*80 pdffile,xsecfile,filemed,filesplit,
buffer,
413 CHARACTER*100 hepmcfile,logfile,filename2
415 LOGICAL pdfexist,splitiexist,xsecexist
417 data maxnscatcen/10000/
426 hepmcfile =
'out.hepmc'
427 filesplit =
'splitint.dat'
429 xsecfile =
'xsecs.dat'
430 filemed =
'medium-params.dat'
453 writescatcen = .
false.
454 writedummies = .
false.
458 if (.not.hadro) shorthepmc = .
true.
464 if (iargc().eq.0)
then
465 write(*,*)
'No parameter file given, '//
466 &
'will run with default settings.'
469 write(*,*)
'Reading parameters from ',
filename
472 read(1,
'(A)', iostat=ios)
buffer
473 if(ios.ne.0) goto 130
475 if (firstchar.eq.
'#') goto 120
479 if(label.eq.
"NEVENT")
then
480 read(
value,*,iostat=ios) nsim
481 elseif(label.eq.
"NJOB")
then
482 read(
value,*,iostat=ios) njob
483 elseif(label.eq.
"LOGFILE")
then
484 read(
value,
'(a)',iostat=ios) logfile
485 elseif(label.eq.
"HEPMCFILE")
then
486 read(
value,
'(a)',iostat=ios) hepmcfile
487 elseif(label.eq.
"SPLITINTFILE")
then
488 read(
value,
'(a)',iostat=ios) filesplit
489 elseif(label.eq.
"PDFFILE")
then
490 read(
value,
'(a)',iostat=ios) pdffile
491 elseif(label.eq.
"XSECFILE")
then
492 read(
value,
'(a)',iostat=ios) xsecfile
493 elseif(label.eq.
"MEDIUMPARAMS")
then
494 read(
value,
'(a)',iostat=ios) filemed
495 elseif(label.eq.
"NF")
then
496 read(
value,*,iostat=ios) nf
497 elseif(label.eq.
"LAMBDAQCD")
then
498 read(
value,*,iostat=ios) lqcd
499 elseif(label.eq.
"Q0")
then
500 read(
value,*,iostat=ios) q0
501 elseif(label.eq.
"PTMIN")
then
502 read(
value,*,iostat=ios) ptmin
503 elseif(label.eq.
"PTMAX")
then
504 read(
value,*,iostat=ios) ptmax
505 elseif(label.eq.
"ETAMAX")
then
507 elseif(label.eq.
"PROCESS")
then
508 read(
value,*,iostat=ios) collider
509 elseif(label.eq.
"ISOCHANNEL")
then
510 read(
value,*,iostat=ios) isochannel
511 elseif(label.eq.
"CHANNEL")
then
513 elseif(label.eq.
"SQRTS")
then
514 read(
value,*,iostat=ios) sqrts
515 elseif(label.eq.
"PDFSET")
then
517 elseif(label.eq.
"NSET")
then
518 read(
value,*,iostat=ios) nset
519 elseif(label.eq.
"MASS")
then
521 elseif(label.eq.
"NPROTON")
then
522 read(
value,*,iostat=ios) nproton
523 elseif(label.eq.
"WEIGHTED")
then
524 read(
value,*,iostat=ios) weighted
525 elseif(label.eq.
"WEXPO")
then
526 read(
value,*,iostat=ios) weightex
527 elseif(label.eq.
"ANGORD")
then
528 read(
value,*,iostat=ios) angord
529 elseif(label.eq.
"KEEPRECOILS")
then
530 read(
value,*,iostat=ios) allhad
531 elseif(label.eq.
"HADRO")
then
532 read(
value,*,iostat=ios) hadro
533 elseif(label.eq.
"HADROTYPE")
then
534 read(
value,*,iostat=ios) hadrotype
535 elseif(label.eq.
"SHORTHEPMC")
then
536 read(
value,*,iostat=ios) shorthepmc
537 elseif(label.eq.
"COMPRESS")
then
539 elseif(label.eq.
"WRITESCATCEN")
then
540 read(
value,*,iostat=ios) writescatcen
541 elseif(label.eq.
"WRITEDUMMIES")
then
542 read(
value,*,iostat=ios) writedummies
544 write(*,*)
'unknown label ',label
550 &
'Unable to open parameter file, will exit the run.'
553 130
close(1,status=
'keep')
557 if (ptmin.lt.3.d0) ptmin = 3.d0
558 if (.not.writescatcen) writedummies = .
false.
560 OPEN(unit=logfid,
file=logfile,status=
'unknown')
568 write(logfid,*)
'parameters of the run:'
569 write(logfid,*)
'NEVENT = ',nsim
570 write(logfid,*)
'NJOB = ',njob
571 write(logfid,*)
'LOGFILE = ',logfile
572 write(logfid,*)
'HEPMCFILE = ',hepmcfile
573 write(logfid,*)
'SPLITINTFILE = ',filesplit
574 write(logfid,*)
'PDFFILE = ',pdffile
575 write(logfid,*)
'XSECFILE = ',xsecfile
576 write(logfid,*)
'MEDIUMPARAMS = ',filemed
577 write(logfid,*)
'NF = ',nf
578 write(logfid,*)
'LAMBDAQCD = ',lqcd
579 write(logfid,*)
'Q0 = ',q0
580 write(logfid,*)
'PTMIN = ',ptmin
581 write(logfid,*)
'PTMAX = ',ptmax
582 write(logfid,*)
'ETAMAX = ',
etamax
583 write(logfid,*)
'PROCESS = ',collider
584 write(logfid,*)
'ISOCHANNEL = ',isochannel
585 write(logfid,*)
'CHANNEL = ',
channel
586 write(logfid,*)
'SQRTS = ',sqrts
587 write(logfid,*)
'PDFSET = ',
pdfset
588 write(logfid,*)
'NSET = ',nset
589 write(logfid,*)
'MASS = ',
mass
590 write(logfid,*)
'NPROTON = ',nproton
591 write(logfid,*)
'WEIGHTED = ',weighted
592 write(logfid,*)
'WEXPO = ',weightex
593 write(logfid,*)
'ANGORD = ',angord
594 write(logfid,*)
'KEEPRECOILS = ',allhad
595 write(logfid,*)
'HADRO = ',hadro
596 write(logfid,*)
'HADROTYPE = ',hadrotype
597 write(logfid,*)
'SHORTHEPMC = ',shorthepmc
598 write(logfid,*)
'COMPRESS = ',
compress
599 write(logfid,*)
'WRITESCATCEN = ',writescatcen
600 write(logfid,*)
'WRITEDUMMIES = ',writedummies
604 if ((collider.ne.
'PPJJ').and.(collider.ne.
'EEJJ')
605 & .and.(collider.ne.
'PPYJ').and.(collider.ne.
'PPYQ')
606 & .and.(collider.ne.
'PPYG')
607 & .and.(collider.ne.
'PPZJ').and.(collider.ne.
'PPZQ')
608 & .and.(collider.ne.
'PPZG').and.(collider.ne.
'PPWJ')
609 & .and.(collider.ne.
'PPWQ').and.(collider.ne.
'PPWG')
610 & .and.(collider.ne.
'PPDY'))
then
611 write(logfid,*)
'Fatal error: colliding system unknown, '//
621 OPEN(unit=hpmcfid,
file=hepmcfile,status=
'unknown')
623 WRITE(hpmcfid,
'(A)')
'HepMC::Version 2.06.05'
624 WRITE(hpmcfid,
'(A)')
'HepMC::IO_GenEvent-START_EVENT_LISTING'
629 eovest=min(1.5*(ptmax+50.)*cosh(
etamax),sqrts/2.)
639 INQUIRE(
file=filesplit,
exist=splitiexist)
641 write(logfid,*)
'read splitting integrals from ',filesplit
642 OPEN(unit=10,
file=filesplit,status=
'old')
643 READ(10,*)qmax,zmmin,npoint
645 READ(10,*) qval(
i),zmval(
i)
649 READ(10,*)splitiggv(
i,
j),splitiqqv(
i,
j),splitiqgv(
i,
j)
652 CLOSE(10,status=
'keep')
654 write(logfid,*)
'have to integrate splitting functions, '//
655 &
'this may take some time'
657 INQUIRE(
file=filesplit,
exist=splitiexist)
658 IF(.NOT.splitiexist)
THEN
659 write(logfid,*)
'write splitting integrals to ',filesplit
660 OPEN(unit=10,
file=filesplit,status=
'new')
661 WRITE(10,*)qmax,zmmin,npoint
663 WRITE(10,*) qval(
i),zmval(
i)
667 WRITE(10,*)splitiggv(
i,
j),splitiqqv(
i,
j),splitiqgv(
i,
j)
670 CLOSE(10,status=
'keep')
677 write(logfid,*)
'read pdfs from ',pdffile
678 OPEN(unit=10,
file=pdffile,status=
'old')
681 READ(10,*)qinqx(
i,
j),ginqx(
i,
j),qingx(
i,
j),gingx(
i,
j)
684 CLOSE(10,status=
'keep')
686 write(logfid,*)
'have to integrate pdfs, this may take some time'
689 IF(.NOT.pdfexist)
THEN
690 write(logfid,*)
'write pdfs to ',pdffile
691 OPEN(unit=10,
file=pdffile,status=
'new')
694 WRITE(10,*)qinqx(
i,
j),ginqx(
i,
j),qingx(
i,
j),gingx(
i,
j)
697 CLOSE(10,status=
'keep')
704 write(logfid,*)
'read cross sections from ',xsecfile
705 OPEN(unit=10,
file=xsecfile,status=
'old')
708 READ(10,*)intq1(
j,jj),intq2(
j,jj),
709 &intg1(
j,jj),intg2(
j,jj)
712 CLOSE(10,status=
'keep')
714 write(logfid,*)
'have to integrate cross sections, '//
715 &
'this may take some time'
718 IF(.NOT.xsecexist)
THEN
719 write(logfid,*)
'write cross sections to ',xsecfile
720 OPEN(unit=10,
file=xsecfile,status=
'new')
723 WRITE(10,*)intq1(
j,jj),intq2(
j,jj),
724 &intg1(
j,jj),intg2(
j,jj)
727 CLOSE(10,status=
'keep')
769 WRITE(snset,
'(i1)') nset
771 WRITE(snset,
'(i2)') nset
773 initstr=
'EPS09LO,'//snset
791 common/
pydat1/mstu(200),paru(200),mstj(200),parj(200)
793 DOUBLE PRECISION paru,parj
794 common/
pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
795 INTEGER mdcy,mdme,kfdp
796 DOUBLE PRECISION brat
797 common/
pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
798 INTEGER msel,mselpd,msub,kfin
799 DOUBLE PRECISION ckin
803 common/
pydatr/mrpy(6),rrpy(100)
805 DOUBLE PRECISION rrpy
807 common/npdf/
mass,nset,eps09,initstr
809 DOUBLE PRECISION mass
816 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
819 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
820 LOGICAL angord,scatrecoil,allhad,
compress
822 common/disc/ndisc,nstrange,ngood,errcount,wdisc,discard
824 INTEGER ndisc,nstrange,ngood,errcount
825 double precision wdisc
827 common/
weight/evweight,sumofweights
828 double precision evweight,sumofweights
830 common/wexpo/weightex
831 DOUBLE PRECISION weightex
836 common/evrecord/nsim,
npart,
offset,hadrotype,sqrts,collider,hadro,
839 double precision sqrts
841 character*2 isochannel
842 logical hadro,shorthepmc
844 common/pythiaparams/ptmin,ptmax,weighted
845 double precision ptmin,ptmax
849 character*2 beam1,beam2
873 IF(collider.EQ.
'PPYQ')
THEN
876 ELSEIF(collider.EQ.
'PPYG')
THEN
880 ELSEIF(collider.EQ.
'PPYJ')
THEN
885 ELSEIF((collider.EQ.
'PPZJ').or.(collider.EQ.
'PPZQ')
886 & .or.(collider.EQ.
'PPZG')
887 & .or.(collider.eq.
'PPDY'))
THEN
889 IF((collider.EQ.
'PPZJ').or.(collider.EQ.
'PPZQ')) msub(30)=1
890 IF((collider.EQ.
'PPZJ').or.(collider.EQ.
'PPZG')) msub(15)=1
891 IF(collider.EQ.
'PPDY') msub(1)=1
909 ELSEIF((collider.EQ.
'PPWJ').or.(collider.EQ.
'PPWQ')
910 & .or.(collider.EQ.
'PPWG'))
THEN
912 IF((collider.EQ.
'PPWJ').or.(collider.EQ.
'PPWQ')) msub(31)=1
913 IF((collider.EQ.
'PPWJ').or.(collider.EQ.
'PPWG')) msub(16)=1
948 IF(weighted)
mstp(142)=1
951 mstu(22)=
max(10,int(5.*nsim/100.))
960 IF(collider.EQ.
'EEJJ')
THEN
962 CALL
pyinit(
'CMS',beam1,beam2,sqrts)
963 ELSEIF((collider.EQ.
'PPJJ').OR.(collider.EQ.
'PPYJ').OR.
964 & (collider.EQ.
'PPYG').OR.(collider.EQ.
'PPYQ'))
THEN
966 CALL
pyinit(
'CMS',beam1,beam2,sqrts)
967 ELSEIF((collider.EQ.
'PPWJ').OR.(collider.EQ.
'PPZJ').or.
968 & (collider.EQ.
'PPWQ').OR.(collider.EQ.
'PPZQ').or.
969 & (collider.EQ.
'PPWG').OR.(collider.EQ.
'PPZG'))
THEN
971 CALL
pyinit(
'CMS',beam1,beam2,sqrts)
972 elseif (collider.eq.
'PPDY')
then
973 CALL
pyinit(
'CMS',beam1,beam2,sqrts)
986 common/hepmcid/hpmcfid,logfid
987 integer hpmcfid,logfid
994 common/
pydat1/mstu(200),paru(200),mstj(200),parj(200)
996 DOUBLE PRECISION paru,parj
997 common/
pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
998 INTEGER mdcy,mdme,kfdp
999 DOUBLE PRECISION brat
1000 common/
pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
1001 INTEGER msel,mselpd,msub,kfin
1002 DOUBLE PRECISION ckin
1006 common/
pydatr/mrpy(6),rrpy(100)
1008 DOUBLE PRECISION rrpy
1010 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
1013 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
1014 LOGICAL angord,scatrecoil,allhad,
compress
1016 common/disc/ndisc,nstrange,ngood,errcount,wdisc,discard
1018 INTEGER ndisc,nstrange,ngood,errcount
1019 double precision wdisc
1021 common/angor/za(23000),zd(23000),thetaa(23000),qqbard(23000)
1022 DOUBLE PRECISION za,zd,thetaa
1025 common/ftimefac/ftfac
1026 DOUBLE PRECISION ftfac
1028 common/time/mv(23000,5)
1031 common/colour/trip(23000),anti(23000),colmax
1032 INTEGER trip,anti,colmax
1034 common/check/nscat,nscateff,nsplit
1035 DOUBLE PRECISION nscat,nscateff,nsplit
1037 common/
weight/evweight,sumofweights
1038 double precision evweight,sumofweights
1040 common/wexpo/weightex
1041 DOUBLE PRECISION weightex
1046 common/jetpoint/
x0,y0
1047 double precision x0,y0
1049 common/evrecord/nsim,
npart,
offset,hadrotype,sqrts,collider,hadro,
1050 &shorthepmc,
channel,isochannel
1052 double precision sqrts
1054 character*2 isochannel
1055 logical hadro,shorthepmc
1057 common/storescatcen/nscatcen,maxnscatcen,scatflav(10000),
1058 &scatcen(10000,5),writescatcen,writedummies
1059 integer nscatcen,maxnscatcen,scatflav
1060 double precision scatcen
1061 logical writescatcen,writedummies
1064 INTEGER nold,pid,ipart,lme1,lme2,
j,
i,lme1orig,lme2orig,llep1,
1066 DOUBLE PRECISION pyr,eni,qmax1,
r,
getmass,
pyp,q1,q2,p21,p22,etot,
1067 &qmax2,pold,en1,en2,
beta(3),enew1,enew2,emax,lambda,x1,x2,x3,
1068 &meweight,psweight,
weight,eps1,eps2,theta1,theta2,z1,z2,
1071 CHARACTER*2 type1,type2
1072 LOGICAL firsttrip,which1,which2,
isdiquark
1073 DATA pi/3.141592653589793d0/
1097 sumofweights=sumofweights+evweight
1098 IF((collider.EQ.
'EEJJ').AND.(abs(
k(8,2)).GT.6))
THEN
1099 wdisc=wdisc+evweight
1107 if (collider.eq.
'PPDY')
then
1115 if((collider.EQ.
'PPZJ').OR.(collider.EQ.
'PPZQ').or.
1116 & (collider.EQ.
'PPZG').or.(collider.EQ.
'PPWJ').or.
1117 & (collider.EQ.
'PPWQ').or.(collider.EQ.
'PPWG'))
THEN
1120 if(abs(
k(7,2)).gt.21)
then
1131 if((collider.EQ.
'PPZJ').OR.(collider.EQ.
'PPZQ').or.
1132 & (collider.EQ.
'PPZG').or.(collider.EQ.
'PPWJ').or.
1133 & (collider.EQ.
'PPWQ').or.(collider.EQ.
'PPWG'))
THEN
1134 if(
k(ipart,3).eq.
offset-1) llep1=ipart
1135 if(
k(ipart,3).eq.
offset) llep2=ipart
1137 IF(
k(ipart,3).EQ.(lme1orig))
THEN
1139 IF(
k(ipart,2).EQ.21)
THEN
1144 ELSEIF(
k(ipart,3).EQ.lme2orig)
THEN
1146 IF(
k(ipart,2).EQ.21)
THEN
1158 IF(
k(ipart,1).EQ.2)
THEN
1159 IF(
k(ipart-1,1).EQ.2)
THEN
1162 trip(ipart)=colmax+1
1163 anti(ipart)=trip(ipart-1)
1165 trip(ipart)=anti(ipart-1)
1166 anti(ipart)=colmax+1
1171 IF(((abs(
k(ipart,2)).LT.10).AND.(
k(ipart,2).GT.0))
1172 & .OR.(
isdiquark(
k(ipart,2)).AND.(
k(ipart,2).LT.0)))
THEN
1173 trip(ipart)=colmax+1
1178 anti(ipart)=colmax+1
1184 IF(
k(ipart,1).EQ.1)
THEN
1188 anti(ipart)=trip(ipart-1)
1190 trip(ipart)=anti(ipart-1)
1195 if (
k(lme1,1).lt.11)
k(lme1,1)=1
1196 if (
k(lme2,1).lt.11)
k(lme2,1)=1
1198 eni=
max(
p(lme1,4),
p(lme2,4))
1200 IF((ipart.NE.lme1).AND.(ipart.NE.lme2).AND.(
k(ipart,1).LT.11))
1202 if (
k(ipart,2).eq.22)
k(ipart,1)=4
1206 if((collider.EQ.
'PPZJ').OR.(collider.EQ.
'PPZQ').or.
1207 & (collider.EQ.
'PPZG').or.(collider.EQ.
'PPWJ').or.
1208 & (collider.EQ.
'PPWQ').or.(collider.EQ.
'PPWG'))
THEN
1209 if (abs(
k(lme1,2)).gt.21)
then
1211 qmax2=sqrt(
pari(18)+
p(lme1,5)**2)
1213 qmax1=sqrt(
pari(18)+
p(lme2,5)**2)
1216 emax=
p(lme1,4)+
p(lme2,4)
1219 ELSEIF(collider.EQ.
'PPJJ'.OR.collider.EQ.
'PPYJ'
1220 & .OR.collider.EQ.
'PPYQ'.OR.collider.EQ.
'PPYG')
THEN
1221 if (
k(lme1,1).eq.4)
then
1226 if (
k(lme2,1).eq.4)
then
1233 emax=
p(lme1,4)+
p(lme2,4)
1239 beta(1)=(
p(lme1,1)+
p(lme2,1))/(
p(lme1,4)+
p(lme2,4))
1240 beta(2)=(
p(lme1,2)+
p(lme2,2))/(
p(lme1,4)+
p(lme2,4))
1241 beta(3)=(
p(lme1,3)+
p(lme2,3))/(
p(lme1,4)+
p(lme2,4))
1244 etot=
p(lme1,4)+
p(lme2,4)
1245 IF(collider.EQ.
'EEJJ')
THEN
1248 emax=
p(lme1,4)+
p(lme2,4)
1257 182
if (abs(
k(lme1,2)).gt.21)
then
1262 if (abs(
k(lme2,2)).gt.21)
then
1267 enew1=etot/2.d0 + (m1**2-m2**2)/(2.*etot)
1268 enew2=etot/2.d0 - (m1**2-m2**2)/(2.*etot)
1269 p21 = (etot/2.d0 + (m1**2-m2**2)/(2.*etot))**2 - m1**2
1270 p22 = (etot/2.d0 - (m1**2-m2**2)/(2.*etot))**2 - m2**2
1272 IF((
pyr(0).GT.
weight).OR.(p21.LT.0.d0).OR.(p22.LT.0.d0)
1273 & .OR.(enew1.LT.0.d0).OR.(enew2.LT.0.d0)
1285 p(lme1,1)=
p(lme1,1)*sqrt(p21)/pold
1286 p(lme1,2)=
p(lme1,2)*sqrt(p21)/pold
1287 p(lme1,3)=
p(lme1,3)*sqrt(p21)/pold
1291 p(lme2,1)=
p(lme2,1)*sqrt(p22)/pold
1292 p(lme2,2)=
p(lme2,2)*sqrt(p22)/pold
1293 p(lme2,3)=
p(lme2,3)*sqrt(p22)/pold
1300 eps1=0.5-0.5*sqrt(1.-q0**2/q1**2)
1301 & *sqrt(1.-q1**2/
p(lme1,4)**2)
1302 IF((z1.LT.eps1).OR.(z1.GT.(1.-eps1)))
THEN
1311 eps2=0.5-0.5*sqrt(1.-q0**2/q2**2)
1312 & *sqrt(1.-q2**2/
p(lme2,4)**2)
1313 IF((z2.LT.eps2).OR.(z2.GT.(1.-eps2)))
THEN
1323 IF(collider.EQ.
'EEJJ')
THEN
1324 beta(1)=(
p(lme1,1)+
p(lme2,1))/(
p(lme1,4)+
p(lme2,4))
1325 beta(2)=(
p(lme1,2)+
p(lme2,2))/(
p(lme1,4)+
p(lme2,4))
1326 beta(3)=(
p(lme1,3)+
p(lme2,3))/(
p(lme1,4)+
p(lme2,4))
1331 x1=z1*(etot**2+q1**2)/etot**2
1332 x2=(etot**2-q1**2)/etot**2
1333 x3=(1.-z1)*(etot**2+q1**2)/etot**2
1334 psweight=(1.-x1)*(1.+(x1/(2.-x2))**2)/x3
1335 & + (1.-x2)*(1.+(x2/(2.-x1))**2)/x3
1336 meweight=x1**2+x2**2
1346 x1=(etot**2-q2**2)/etot**2
1347 x2=z2*(etot**2+q2**2)/etot**2
1348 x3=(1.-z2)*(etot**2+q2**2)/etot**2
1349 psweight=(1.-x1)*(1.+(x1/(2.-x2))**2)/x3
1350 & + (1.-x2)*(1.+(x2/(2.-x1))**2)/x3
1351 meweight=x1**2+x2**2
1358 186 enew1=etot/2.d0 + (q1**2-q2**2)/(2.*etot)
1359 enew2=etot/2.d0 - (q1**2-q2**2)/(2.*etot)
1360 p21 = (etot/2.d0 + (q1**2-q2**2)/(2.*etot))**2 - q1**2
1361 p22 = (etot/2.d0 - (q1**2-q2**2)/(2.*etot))**2 - q2**2
1363 p(lme1,1)=
p(lme1,1)*sqrt(p21)/pold
1364 p(lme1,2)=
p(lme1,2)*sqrt(p21)/pold
1365 p(lme1,3)=
p(lme1,3)*sqrt(p21)/pold
1369 p(lme2,1)=
p(lme2,1)*sqrt(p22)/pold
1370 p(lme2,2)=
p(lme2,2)*sqrt(p22)/pold
1371 p(lme2,3)=
p(lme2,3)*sqrt(p22)/pold
1378 eps1=0.5-0.5*sqrt(1.-q0**2/q1**2)
1379 & *sqrt(1.-q1**2/
p(lme1,4)**2)
1380 IF((z1.LT.eps1).OR.(z1.GT.(1.-eps1)))
THEN
1389 eps2=0.5-0.5*sqrt(1.-q0**2/q2**2)
1390 & *sqrt(1.-q2**2/
p(lme2,4)**2)
1391 IF((z2.LT.eps2).OR.(z2.GT.(1.-eps2)))
THEN
1402 if((collider.EQ.
'PPZJ').OR.(collider.EQ.
'PPZQ').or.
1403 & (collider.EQ.
'PPZG').or.(collider.EQ.
'PPWJ').or.
1404 & (collider.EQ.
'PPWQ').or.(collider.EQ.
'PPWG'))
THEN
1410 if (abs(
k(lme1,2)).gt.21)
then
1411 beta(1)=
p(lme1,1)/
p(lme1,4)
1412 beta(2)=
p(lme1,2)/
p(lme1,4)
1413 beta(3)=
p(lme1,3)/
p(lme1,4)
1415 beta(1)=
p(lme2,1)/
p(lme2,4)
1416 beta(2)=
p(lme2,2)/
p(lme2,4)
1417 beta(3)=
p(lme2,3)/
p(lme2,4)
1426 thetaa(lme1)=
p(lme1,5)/(sqrt(z1*(1.-z1))*
p(lme1,4))
1427 thetaa(lme2)=
p(lme2,5)/(sqrt(z2*(1.-z2))*
p(lme2,4))
1437 IF(
p(lme1,5).GT.0.d0)
THEN
1438 lambda=1.d0/(ftfac*
p(lme1,4)*0.2/q1**2)
1439 mv(lme1,5)=-
log(1.d0-
pyr(0))/lambda
1448 IF(
p(lme2,5).GT.0.d0)
THEN
1449 lambda=1.d0/(ftfac*
p(lme2,4)*0.2/q2**2)
1450 mv(lme2,5)=-
log(1.d0-
pyr(0))/lambda
1459 wdisc=wdisc+evweight
1461 write(logfid,*)
'discard event',
j
1467 IF(
k(
i,1).EQ.3)
k(
i,1)=22
1473 write(logfid,*)
'discard event',
j
1474 wdisc=wdisc+evweight
1480 IF(mstu(30).NE.errcount)
THEN
1481 write(logfid,*)
'PYTHIA discards event',
j,
1482 &
' (error number',mstu(30),
')'
1484 wdisc=wdisc+evweight
1491 IF(mstu(30).NE.errcount)
THEN
1498 102
IF(nsim.GT.100)
THEN
1499 IF(mod(
j,nsim/100).EQ.0)
THEN
1500 write(logfid,*)
'done with event number ',
j
1503 write(logfid,*)
'done with event number ',
j
1516 common/hepmcid/hpmcfid,logfid
1517 integer hpmcfid,logfid
1521 ELSEIF(which.EQ.1)
THEN
1524 WRITE(logfid,*)
'error: unknown hadronisation type in MAKESTRINGS'
1535 common/hepmcid/hpmcfid,logfid
1536 integer hpmcfid,logfid
1540 DOUBLE PRECISION p,
v
1542 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
1545 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
1546 LOGICAL angord,scatrecoil,allhad,
compress
1548 common/colour/trip(23000),anti(23000),colmax
1549 INTEGER trip,anti,colmax
1551 common/disc/ndisc,nstrange,ngood,errcount,wdisc,discard
1553 INTEGER ndisc,nstrange,ngood,errcount
1554 double precision wdisc
1556 INTEGER nold,
i,
j,lquark,lmatch,lloose,nold1
1557 DOUBLE PRECISION eaddend,
pyr,
dir
1566 IF((
k(
i,1).EQ.4).AND.(trip(
i).EQ.0).AND.(anti(
i).EQ.0))
THEN
1570 write(logfid,*)
'event too long for event record'
1596 IF((
k(
i,1).EQ.11).OR.(
k(
i,1).EQ.12).OR.(
k(
i,1).EQ.13)
1597 & .OR.(
k(
i,1).EQ.14))
k(
i,1)=17
1598 IF(((
k(
i,1).EQ.1).OR.(
k(
i,1).EQ.3).OR.(
k(
i,1).EQ.4).OR.
1609 write(logfid,*)
'event too long for event record'
1626 trip(
n)=trip(lquark)
1627 anti(
n)=anti(lquark)
1632 IF(((
k(
i,1).EQ.1).OR.(
k(
i,1).EQ.3).OR.(
k(
i,1).EQ.4)
1633 & .OR.(
k(
i,1).EQ.5))
1634 & .AND.(((trip(
i).EQ.anti(
n)).AND.(trip(
i).NE.0))
1635 & .OR.((anti(
i).EQ.trip(
n)).AND.(anti(
i).NE.0))))
THEN
1638 write(logfid,*)
'event too long for event record'
1656 IF(
k(
i,2).EQ.21)
THEN
1666 write(logfid,*)
'Error in MAKESTRINGS_VAC: failed to reconstruct '//
1667 &
'colour singlet system, will discard event'
1676 IF(((
k(
i,1).EQ.1).OR.(
k(
i,1).EQ.3).OR.(
k(
i,1).EQ.4)
1677 & .OR.(
k(
i,1).EQ.5)))
THEN
1679 IF(((
k(
i,1).EQ.1).OR.(
k(
i,1).EQ.3).OR.(
k(
i,1).EQ.4)
1680 & .OR.(
k(
i,1).EQ.5)))
THEN
1681 IF(anti(
i).EQ.trip(
j)) goto 45
1691 write(logfid,*)
'Error in MAKESTRINGS_VAC: failed to reconstruct '//
1692 &
'colour singlet system, will discard event'
1698 write(logfid,*)
'event too long for event record'
1715 trip(
n)=trip(lloose)
1716 anti(
n)=anti(lloose)
1721 IF(((
k(
i,1).EQ.1).OR.(
k(
i,1).EQ.3).OR.(
k(
i,1).EQ.4)
1722 & .OR.(
k(
i,1).EQ.5))
1723 & .AND.(anti(
i).EQ.trip(
n)))
THEN
1726 write(logfid,*)
'event too long for event record'
1749 write(logfid,*)
'Error in MAKESTRINGS_VAC: failed to reconstruct '//
1750 &
'colour singlet system, will discard event'
1767 DOUBLE PRECISION p,
v
1769 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
1772 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
1773 LOGICAL angord,scatrecoil,allhad,
compress
1775 common/colour/trip(23000),anti(23000),colmax
1776 INTEGER trip,anti,colmax
1778 INTEGER nold,
i,
j,lmax,lmin,lend,nold1
1784 common/hepmcid/hpmcfid,logfid
1785 integer hpmcfid,logfid
1787 common/disc/ndisc,nstrange,ngood,errcount,wdisc,discard
1789 INTEGER ndisc,nstrange,ngood,errcount
1790 double precision wdisc
1798 IF((
k(
i,1).EQ.4).AND.(trip(
i).EQ.0).AND.(anti(
i).EQ.0))
THEN
1802 write(logfid,*)
'event too long for event record'
1828 IF((
k(
i,1).EQ.11).OR.(
k(
i,1).EQ.12).OR.(
k(
i,1).EQ.13)
1829 & .OR.(
k(
i,1).EQ.14))
k(
i,1)=17
1830 if (abs(
pyp(
i,17)).gt.4.d0)
k(
i,1)=17
1831 IF(((
k(
i,1).EQ.1).OR.(
k(
i,1).EQ.3).OR.(
k(
i,1).EQ.4)
1832 & .OR.(
k(
i,1).EQ.5)).AND.(
p(
i,4).GT.emax))
THEN
1838 IF(lmax.EQ.0) goto 50
1840 IF(
k(lmax,2).EQ.21)
THEN
1845 write(logfid,*)
'event too long for event record'
1849 IF((
n-2).GT.nold)
THEN
1851 k(
n+nold-
j,1)=
k(
n+nold-
j-2,1)
1852 k(
n+nold-
j,2)=
k(
n+nold-
j-2,2)
1853 IF(
k(
n+nold-
j-2,3).GT.nold)
THEN
1854 k(
n+nold-
j,3)=
k(
n+nold-
j-2,3)+2
1856 k(
n+nold-
j,3)=
k(
n+nold-
j-2,3)
1860 p(
n+nold-
j,1)=
p(
n+nold-
j-2,1)
1861 p(
n+nold-
j,2)=
p(
n+nold-
j-2,2)
1862 p(
n+nold-
j,3)=
p(
n+nold-
j-2,3)
1863 p(
n+nold-
j,4)=
p(
n+nold-
j-2,4)
1864 p(
n+nold-
j,5)=
p(
n+nold-
j-2,5)
1865 k(
k(
n+nold-
j-2,3),4)=
k(
k(
n+nold-
j-2,3),4)+2
1866 k(
k(
n+nold-
j-2,3),5)=
k(
k(
n+nold-
j-2,3),5)+2
1884 p(nold-1,1)=(1.-
z)*
p(lmax,1)
1885 p(nold-1,2)=(1.-
z)*
p(lmax,2)
1886 p(nold-1,3)=(1.-
z)*
p(lmax,3)
1887 p(nold-1,4)=(1.-
z)*
p(lmax,4)
1888 p(nold-1,5)=
p(lmax,5)
1893 p(nold,1)=
z*
p(lmax,1)
1894 p(nold,2)=
z*
p(lmax,2)
1895 p(nold,3)=
z*
p(lmax,3)
1896 p(nold,4)=
z*
p(lmax,4)
1905 write(logfid,*)
'event too long for event record'
1927 IF(((
k(
i,1).EQ.1).OR.(
k(
i,1).EQ.3).OR.(
k(
i,1)
1928 & .EQ.4).OR.(
k(
i,1).EQ.5))
1929 & .AND.((
k(
i,2).EQ.21).OR.((
k(
i,2)*
k(lend,2).LT.0.d0).AND.
1930 & (
k(
i,3).NE.
k(lend,3))))
1931 & .AND.(
p(
i,1)*
p(lend,1).GT.0.d0))
THEN
1932 minv=
p(
i,4)*
p(lmax,4)-
p(
i,1)*
p(lmax,1)-
p(
i,2)*
p(lmax,2)
1934 IF((minv.LT.mmin).AND.(minv.GT.0.d0).AND.(minv.LT.mcut))
THEN
1944 write(logfid,*)
'event too long for event record'
1955 IF(
pyr(0).LT.0.5)
THEN
1968 write(logfid,*)
'event too long for event record'
1984 IF(
k(lmin,2).EQ.21)
THEN
2006 DOUBLE PRECISION p,
v
2008 INTEGER nfirst,nlast,
i,
j
2011 DO 21
i=1,nlast-nfirst
2031 DOUBLE PRECISION p,
v
2033 common/time/mv(23000,5)
2036 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
2039 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
2040 LOGICAL angord,scatrecoil,allhad,
compress
2042 common/disc/ndisc,nstrange,ngood,errcount,wdisc,discard
2044 INTEGER ndisc,nstrange,ngood,errcount
2045 double precision wdisc
2056 IF((
k(
i,1).EQ.1).OR.(
k(
i,1).EQ.2))
THEN
2074 DOUBLE PRECISION p,
v
2076 common/time/mv(23000,5)
2079 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
2082 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
2083 LOGICAL angord,scatrecoil,allhad,
compress
2085 common/disc/ndisc,nstrange,ngood,errcount,wdisc,discard
2087 INTEGER ndisc,nstrange,ngood,errcount
2088 double precision wdisc
2090 common/angor/za(23000),zd(23000),thetaa(23000),qqbard(23000)
2091 DOUBLE PRECISION za,zd,thetaa
2094 common/check/nscat,nscateff,nsplit
2095 DOUBLE PRECISION nscat,nscateff,nsplit
2097 common/
coherent/nstart,nend,allqs(10000,6),scatcentres(10000,10),
2100 DOUBLE PRECISION allqs,scatcentres,qsumvec,qsum2
2102 common/
weight/evweight,sumofweights
2103 double precision evweight,sumofweights
2105 common/hepmcid/hpmcfid,logfid
2106 integer hpmcfid,logfid
2108 common/storescatcen/nscatcen,maxnscatcen,scatflav(10000),
2109 & scatcen(10000,5),writescatcen,writedummies
2110 integer nscatcen,maxnscatcen,scatflav
2111 double precision scatcen
2112 logical writescatcen,writedummies
2114 INTEGER l,
line,nold,typi,lineold,lkine,nendold,nscatcenold
2115 DOUBLE PRECISION theta,
phi,
pyp,formtime,starttime,tleft,
2117 LOGICAL overq0,qqbardec
2124 starttime=mv(
line,4)
2139 20
IF(discard)
RETURN
2140 IF(((
k(
line,1).EQ.1).AND.(
p(
line,5).GT.0.d0))
2141 & .OR.((
k(
line,1).EQ.2).AND.(zd(
line).gt.0.d0)))
THEN
2145 formtime=min(mv(
line,5),ltime)
2152 tleft=formtime-starttime
2153 IF(
k(
line,2).EQ.21)
THEN
2160 IF(tleft.LE.1.d-10)
THEN
2166 nscatcenold=nscatcen
2176 CALL
pyrobo(
n-1,
n,theta,0d0,0d0,0d0,0d0)
2179 mv(
n-1,1)=mv(
line,1)
2181 mv(
n-1,2)=mv(
line,2)
2183 mv(
n-1,3)=mv(
line,3)
2228 qqbardec=qqbard(
line)
2229 nscatcenold=nscatcen
2230 25
IF(
x.LT.1.d0)
THEN
2236 IF(
k(
line,2).EQ.21)
THEN
2237 newmass=
getmass(0.d0,scalefacm*sqrt(-qsum2),-1.d0,
p(
line,4),
2238 &
'GC',sqrt(-qsum2),.
false.,zdec,qqbardec)
2239 IF(zdec.GT.0.d0)
THEN
2240 thetaa(
line)=newmass/(sqrt(zdec*(1.-zdec))*
p(
line,4))
2245 qqbard(
line)=qqbardec
2247 newmass=
getmass(0.d0,scalefacm*sqrt(-qsum2),-1.d0,
p(
line,4),
2248 &
'QQ',sqrt(-qsum2),.
false.,zdec,qqbardec)
2249 IF(zdec.GT.0.d0)
THEN
2250 thetaa(
line)=newmass/(sqrt(zdec*(1.-zdec))*
p(
line,4))
2255 qqbard(
line)=qqbardec
2258 qqbardec=qqbard(
line)
2263 qsumvec(1)=allqs(nend,2)
2264 qsumvec(2)=allqs(nend,3)
2265 qsumvec(3)=allqs(nend,4)
2266 qsumvec(4)=allqs(nend,5)
2267 IF(-allqs(nend,1).GT.q0**2/scalefacm**2)
THEN
2272 tleft = starttime+tsum+tleft-allqs(1,6)
2273 tsum = allqs(1,6)-starttime
2280 & newmass,overq0,zdec,qqbardec)
2281 IF(newmass.GT.(
p(
line,5)*(1.d0+1.d-6)))
THEN
2286 qqbardec=qqbard(
line)
2297 CALL
dokinematics(lkine,lineold,nstart,nend,newmass,retrysplit,
2298 & starttime+tsum,
x,zdec,qqbardec)
2300 tleft = starttime+tsum+tleft-allqs(1,6)
2301 tsum = allqs(1,6)-starttime
2305 qsumvec(1)=allqs(nend,2)
2306 qsumvec(2)=allqs(nend,3)
2307 qsumvec(3)=allqs(nend,4)
2308 qsumvec(4)=allqs(nend,5)
2310 IF(-allqs(nend,1).GT.q0**2/scalefacm**2)
THEN
2319 nscatcen=nscatcenold
2320 nsplit=nsplit-evweight
2324 starttime=starttime+tsum
2329 starttime=starttime+tsum
2333 starttime=starttime+tsum
2339 21
IF(((
k(
line,1).EQ.1).AND.(
p(
line,5).GT.0.d0))
2340 & .OR.((
k(
line,1).EQ.2).AND.(zd(
line).gt.0.d0))
2341 & .OR.(starttime.LT.ltime))
THEN
2356 common/hepmcid/hpmcfid,logfid
2357 integer hpmcfid,logfid
2361 DOUBLE PRECISION p,
v
2363 common/time/mv(23000,5)
2366 common/ftimefac/ftfac
2367 DOUBLE PRECISION ftfac
2369 common/colour/trip(23000),anti(23000),colmax
2370 INTEGER trip,anti,colmax
2372 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
2375 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
2376 LOGICAL angord,scatrecoil,allhad,
compress
2378 common/disc/ndisc,nstrange,ngood,errcount,wdisc,discard
2380 INTEGER ndisc,nstrange,ngood,errcount
2381 double precision wdisc
2383 common/angor/za(23000),zd(23000),thetaa(23000),qqbard(23000)
2384 DOUBLE PRECISION za,zd,thetaa
2387 common/check/nscat,nscateff,nsplit
2388 DOUBLE PRECISION nscat,nscateff,nsplit
2390 common/
weight/evweight,sumofweights
2391 double precision evweight,sumofweights
2396 &
getmass,
pz,
eps,qh,
z,
r,lambda,
weight,zdecb,zdecc,xdec(3),theta,
2398 LOGICAL quark,qqbar,qqbardecb,qqbardecc
2400 DATA pi/3.141592653589793d0/
2402 IF((
n+2).GT.22990)
THEN
2403 write(logfid,*)
'event too long for event record'
2408 xdec(1)=mv(l,1)+(mv(l,5)-mv(l,4))*
p(l,1)/
p(l,4)
2409 xdec(2)=mv(l,2)+(mv(l,5)-mv(l,4))*
p(l,2)/
p(l,4)
2410 xdec(3)=mv(l,3)+(mv(l,5)-mv(l,4))*
p(l,3)/
p(l,4)
2411 IF(
gettemp(xdec(1),xdec(2),xdec(3),mv(l,5)).GT.0.d0)
THEN
2418 IF((
p(l,5).EQ.0d0).OR.(
k(l,1).EQ.11).OR.(
k(l,1).EQ.12)
2419 & .OR.(
k(l,1).EQ.13).OR.(
k(l,1).EQ.14).OR.(
k(l,1).EQ.3)
2420 & .or.(zd(l).lt.0.d0)) goto 31
2422 IF(
k(l,2).EQ.21)
THEN
2431 IF(quark.OR.qqbar)
THEN
2434 IF(
pyr(0).LT.0.5)
THEN
2442 write(logfid,*)
'makesplitting: z=0',l
2447 36
IF(angord.AND.(za(l).NE.1.d0))
THEN
2449 qh=4.*
p(l,5)**2*(1.-za(l))/(za(l)*
p(
k(l,3),5)**2)
2451 write(logfid,*)l,
': reject event: angular ordering
2452 & conflict in medium'
2457 eps=0.5-0.5*sqrt(1.-qh)
2472 bmax1=min(
p(l,5),
z*
p(l,4))
2473 cmax1=min(
p(l,5),(1.-
z)*
p(l,4))
2475 30
IF(quark.OR.qqbar)
THEN
2476 mb=
getmass(0.d0,bmax1,theta,
z*
p(l,4),
'QQ',
2477 & bmax1,.
false.,zdecb,qqbardecb)
2479 mb=
getmass(0.d0,bmax1,theta,
z*
p(l,4),
'GC',
2480 & bmax1,.
false.,zdecb,qqbardecb)
2483 IF(quark.OR.(.NOT.qqbar))
THEN
2484 mc=
getmass(0.d0,cmax1,theta,(1.-
z)*
p(l,4),
'GC',
2485 & cmax1,.
false.,zdecc,qqbardecc)
2487 mc=
getmass(0.d0,cmax1,theta,(1.-
z)*
p(l,4),
'QQ',
2488 & cmax1,.
false.,zdecc,qqbardecc)
2491 182
pz=(2.*
z*
p(l,4)**2-
p(l,5)**2-mb**2+mc**2)/(2.*
p(l,3))
2492 pts=
z**2*(
p(l,4)**2)-
pz**2-mb**2
2496 IF((mb.EQ.0.d0).AND.(mc.EQ.0.d0).AND.(pts.LT.0.d0)) goto 36
2499 & .OR.((mb+mc).GT.
p(l,5)))
THEN
2501 IF(quark.OR.qqbar)
THEN
2503 & bmax1,.
false.,zdecb,qqbardecb)
2506 & bmax1,.
false.,zdecb,qqbardecb)
2509 IF(quark.OR.(.NOT.qqbar))
THEN
2510 mc=
getmass(0.d0,mc,theta,(1.-
z)*
p(l,4),
'GC',
2511 & cmax1,.
false.,zdecc,qqbardecc)
2513 mc=
getmass(0.d0,mc,theta,(1.-
z)*
p(l,4),
'QQ',
2514 & cmax1,.
false.,zdecc,qqbardecc)
2528 IF((
k(l,2).GT.0).AND.(
dir.GE.0))
THEN
2540 p(
n-1,4)=(1-
z)*
p(l,4)
2543 IF(zdecc.GT.0.d0)
THEN
2544 thetaa(
n-1)=
p(
n-1,5)/(sqrt(zdecc*(1.-zdecc))*
p(
n-1,4))
2549 qqbard(
n-1)=qqbardecc
2583 IF(zdecb.GT.0.d0)
THEN
2584 thetaa(
n)=
p(
n,5)/(sqrt(zdecb*(1.-zdecb))*
p(
n,4))
2592 p(
n,1)=sqrt(pts)*cos(phiq)
2593 p(
n,2)=sqrt(pts)*sin(phiq)
2595 p(
n-1,1)=
p(l,1)-
p(
n,1)
2596 p(
n-1,2)=
p(l,2)-
p(
n,2)
2597 p(
n-1,3)=
p(l,3)-
p(
n,3)
2599 IF(
p(
n-1,5).GT.0.d0)
THEN
2600 lambda=1.d0/(ftfac*
p(
n-1,4)*0.2/
p(
n-1,5)**2)
2601 mv(
n-1,5)=mv(l,5)-
log(1.d0-
pyr(0))/lambda
2606 IF(
p(
n,5).GT.0.d0)
THEN
2607 lambda=1.d0/(ftfac*
p(
n,4)*0.2/
p(
n,5)**2)
2608 mv(
n,5)=mv(l,5)-
log(1.d0-
pyr(0))/lambda
2620 nsplit=nsplit+evweight
2631 common/hepmcid/hpmcfid,logfid
2632 integer hpmcfid,logfid
2636 DOUBLE PRECISION p,
v
2638 common/time/mv(23000,5)
2641 common/ftimefac/ftfac
2642 DOUBLE PRECISION ftfac
2644 common/colour/trip(23000),anti(23000),colmax
2645 INTEGER trip,anti,colmax
2647 common/angor/za(23000),zd(23000),thetaa(23000),qqbard(23000)
2648 DOUBLE PRECISION za,zd,thetaa
2651 common/disc/ndisc,nstrange,ngood,errcount,wdisc,discard
2653 INTEGER ndisc,nstrange,ngood,errcount
2654 double precision wdisc
2656 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
2659 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
2660 LOGICAL angord,scatrecoil,allhad,
compress
2662 common/check/nscat,nscateff,nsplit
2663 DOUBLE PRECISION nscat,nscateff,nsplit
2665 common/
weight/evweight,sumofweights
2666 double precision evweight,sumofweights
2669 INTEGER l,typi,nold,
dir
2670 DOUBLE PRECISION x,virt,mb2,mc2,
getmass,
pz,
kt2,theta,
phi,
pi,
2671 &phiq,
pyp,
pyr,
r,time,tsum,taurad,lambda,zdec
2673 CHARACTER*2 typ2,typc
2675 DATA pi/3.141592653589793d0/
2677 IF((
n+2).GT.22990)
THEN
2678 write(logfid,*)
'event too long for event record'
2683 IF(
k(l,2).EQ.21)
THEN
2702 IF(typ2.EQ.
'GG')
THEN
2703 IF(
pyr(0).LT.0.5)
THEN
2714 mc2=
getmass(0.d0,scalefacm*sqrt(-tsum),-1.d0,
2715 & (1.-
x)*
p(l,4),typc,(1.-
x)*
p(l,4),
2716 & .
false.,zdec,qqbardec)**2
2723 CALL
pyrobo(l,l,-theta,0d0,0d0,0d0,0d0)
2724 pz=(2*
x*
p(l,4)**2-
p(l,5)**2-mb2+mc2)/(2*
p(l,3))
2725 kt2=
x**2*(
p(l,4)**2)-
pz**2-mb2
2728 pz=(2*
x*
p(l,4)**2-
p(l,5)**2-mb2+mc2)/(2*
p(l,3))
2729 kt2=
x**2*(
p(l,4)**2)-
pz**2-mb2
2731 CALL
pyrobo(l,l,theta,0d0,0d0,0d0,0d0)
2740 IF(typ2.EQ.
'QG')
THEN
2742 IF(
k(
n-1,2).GT.0)
THEN
2749 ELSEIF(typ2.EQ.
'GQ')
THEN
2751 IF(
k(
n-1,2).GT.0)
THEN
2761 IF((
k(l,2).GT.0).AND.(
dir.GE.0))
THEN
2773 p(
n-1,4)=(1.-
x)*
p(l,4)
2778 IF(typ2.EQ.
'QG')
THEN
2787 ELSEIF(typi.NE.21)
THEN
2798 IF(
k(
n-1,2).EQ.21)
THEN
2806 ELSEIF(
k(
n-1,2).GT.0)
THEN
2826 p(
n,1)=sqrt(
kt2)*cos(phiq)
2827 p(
n,2)=sqrt(
kt2)*sin(phiq)
2829 p(
n-1,1)=
p(l,1)-
p(
n,1)
2830 p(
n-1,2)=
p(l,2)-
p(
n,2)
2831 p(
n-1,3)=
p(l,3)-
p(
n,3)
2834 IF(
p(
n-1,5).GT.0.d0)
THEN
2835 lambda=1.d0/(ftfac*
p(
n-1,4)*0.2/
p(
n-1,5)**2)
2836 mv(
n-1,5)=mv(l,5)-
log(1.d0-
pyr(0))/lambda
2841 IF(
p(
n,5).GT.0.d0)
THEN
2849 qqbard(
n-1)=qqbardec
2862 nsplit=nsplit+evweight
2863 CALL
pyrobo(l,l,theta,0d0,0d0,0d0,0d0)
2864 CALL
pyrobo(
n-1,
n,theta,0d0,0d0,0d0,0d0)
2869 mv(
n-1,1)=mv(l,1)+(mv(
n-1,4)-mv(l,4))*
p(l,1)/
max(
pyp(l,8),
p(l,4))
2870 mv(
n-1,2)=mv(l,2)+(mv(
n-1,4)-mv(l,4))*
p(l,2)/
max(
pyp(l,8),
p(l,4))
2871 mv(
n-1,3)=mv(l,3)+(mv(
n-1,4)-mv(l,4))*
p(l,3)/
max(
pyp(l,8),
p(l,4))
2872 mv(
n, 1)=mv(l,1)+(mv(
n, 4)-mv(l,4))*
p(l,1)/
max(
pyp(l,8),
p(l,4))
2873 mv(
n, 2)=mv(l,2)+(mv(
n, 4)-mv(l,4))*
p(l,2)/
max(
pyp(l,8),
p(l,4))
2874 mv(
n, 3)=mv(l,3)+(mv(
n, 4)-mv(l,4))*
p(l,3)/
max(
pyp(l,8),
p(l,4))
2888 DOUBLE PRECISION p,
v
2890 common/time/mv(23000,5)
2893 common/ftimefac/ftfac
2894 DOUBLE PRECISION ftfac
2896 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
2899 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
2900 LOGICAL angord,scatrecoil,allhad,
compress
2902 common/disc/ndisc,nstrange,ngood,errcount,wdisc,discard
2904 INTEGER ndisc,nstrange,ngood,errcount
2905 double precision wdisc
2907 common/
coherent/nstart,nend,allqs(10000,6),scatcentres(10000,10),
2910 DOUBLE PRECISION allqs,scatcentres,qsumvec,qsum2
2912 common/hepmcid/hpmcfid,logfid
2913 integer hpmcfid,logfid
2915 INTEGER l,typi,
counter,countmax,count2
2917 &
weight,
low,fmax,
getpdf,sigmatot,
getsscat,pfchange,
pi,tnow,tleft,
2918 &
xmax,
pqq,
pqg,
pgq,
pgg,
alphas,tstart,tsum,
q,qold,q2old,
getnewmass,
2919 &
generatez,tmax,tmaxnew,
dt,xsc,ysc,zsc,tsc,ms1,md1,
getms,
getmd,
2922 LOGICAL fchange,norad,overq0,noscat,
getdeltat,retrysplit,
2926 DATA pi/3.141592653589793d0/
2927 DATA countmax/10000/
2931 xsc=mv(l,1)+(tstart-mv(l,4))*
p(l,1)/
p(l,4)
2932 ysc=mv(l,2)+(tstart-mv(l,4))*
p(l,2)/
p(l,4)
2933 zsc=mv(l,3)+(tstart-mv(l,4))*
p(l,3)/
p(l,4)
2935 md1=
getmd(xsc,ysc,zsc,tsc)
2936 ms1=
getms(xsc,ysc,zsc,tsc)
2938 IF(md1.LE.1.d-4.OR.ms1.LE.1.d-4)
THEN
2939 write(logfid,*)
'problem!',
gettemp(xsc,ysc,zsc,tsc),
2945 IF(noscat.AND.(.NOT.retrysplit)) goto 116
2949 IF((
pyr(0).LT.pnorad).OR.(
p(l,4).LT.1.001*q0))
THEN
2956 IF(
k(l,2).EQ.21)
THEN
2959 sigmatot=
getsscat(
p(l,4),
p(l,1),
p(l,2),
p(l,3),
p(l,5),
2960 & q0,
'G',
'C',xsc,ysc,zsc,tsc,0)
2961 IF((sigmatot.EQ.0.d0).OR.(pnorad.EQ.1.d0))
THEN
2964 pfchange=
getsscat(
p(l,4),
p(l,1),
p(l,2),
p(l,3),
p(l,5),
2965 & q0,
'G',
'Q',xsc,ysc,zsc,tsc,0)
2968 sigmatot=
getsscat(
p(l,4),
p(l,1),
p(l,2),
p(l,3),
p(l,5),
2969 & 0.d0,
'G',
'C',xsc,ysc,zsc,tsc,0)
2973 sigmatot=
getsscat(
p(l,4),
p(l,1),
p(l,2),
p(l,3),
p(l,5),
2974 & q0,
'Q',
'C',xsc,ysc,zsc,tsc,0)
2975 IF((sigmatot.EQ.0.d0).OR.(pnorad.EQ.1.d0))
THEN
2978 pfchange=
getsscat(
p(l,4),
p(l,1),
p(l,2),
p(l,3),
p(l,5),
2979 & q0,
'Q',
'G',xsc,ysc,zsc,tsc,0)
2982 sigmatot=
getsscat(
p(l,4),
p(l,1),
p(l,2),
p(l,3),
p(l,5),
2983 & 0.d0,
'Q',
'C',xsc,ysc,zsc,tsc,0)
2985 IF((pfchange.LT.-1.d-4).OR.(pfchange.GT.1.d0+1.d-4))
THEN
2986 write(logfid,*)
'error: flavour change probability=',
2987 & pfchange,
'for ',typ
2989 IF(
pyr(0).LT.pfchange)
THEN
2994 IF (norad) fchange=.
false.
2998 typi=int(
sign(2.d0,
pyr(0)-0.5))
3009 low=q0**2/scalefacm**2
3010 tmax=4.*(
p(l,4)**2-
p(l,5)**2)
3011 xmax=1.-q0**2/(scalefacm**2*4.*tmax)
3013 IF(sigmatot.EQ.0.d0) goto 116
3035 tmaxnew=(
x*
p(l,4))**2
3059 118 CALL
getqvec(l,nend,tnow-mv(l,4),
x)
3060 IF(-qsum2.GT.
p(l,4)**2)
THEN
3066 IF(count2.LT.100)
THEN
3078 IF(-allqs(nend,1).GT.
low) overq0=.
true.
3080 IF(overq0.AND.(.NOT.norad))
THEN
3081 q=
getnewmass(l,scalefacm**2*qsum2,scalefacm**2*q2old,0.d0,
3082 & .
true.,
x,zdum,qqbardum)
3088 111
IF((
q.EQ.0.d0).OR.(
q.EQ.
p(l,5)))
THEN
3091 tauest=ftfac*(1.-
phi)*0.2*
x*
p(l,4)/
q**2
3094 tauest=-
log(1.d0-rtau)/lambda
3097 noscat=.NOT.
getdeltat(l,tnow,min(tleft,tauest),deltal)
3098 IF((.NOT.noscat).AND.(.NOT.norad))
THEN
3109 IF((
q.NE.0.d0).AND.(
q.NE.
p(l,5)))
THEN
3115 qsumvecold(1)=qsumvec(1)
3116 qsumvecold(2)=qsumvec(2)
3117 qsumvecold(3)=qsumvec(3)
3118 qsumvecold(4)=qsumvec(4)
3120 119 CALL
getqvec(l,nend,tnow-mv(l,4),
x)
3121 IF(-qsum2.GT.
p(l,4)**2)
THEN
3122 qsumvec(1)=qsumvecold(1)
3123 qsumvec(2)=qsumvecold(2)
3124 qsumvec(3)=qsumvecold(3)
3125 qsumvec(4)=qsumvecold(4)
3127 IF(count2.LT.100)
THEN
3140 & .OR.(-allqs(nend,1).GT.
low)) overq0=.
true.
3143 IF(overq0.AND.(.NOT.norad))
THEN
3144 q=
getnewmass(l,scalefacm**2*qsum2,scalefacm**2*q2old,0.d0,
3145 & .
true.,
x,zdum,qqbardum)
3153 114 tmaxnew=
x**2*
p(l,4)**2
3158 ELSEIF((-qsum2.LT.
low).OR.(
q.EQ.0.d0))
THEN
3160 ELSEIF(-qsum2.GT.
p(l,4)**2)
THEN
3164 fmax=2.*
log(-scalefacm**2*qsum2/q0**2)
3166 IF(qsum2.EQ.0.d0)
THEN
3173 write(logfid,*)
'x,sqrt(qsum^2),getpdf,fmax:',
x,
3174 & sqrt(-qsum2),
getpdf(
x,scalefacm*sqrt(-qsum2),
'QG'),
'qg',
3180 write(logfid,*)
'x,sqrt(qsum^2),getpdf,fmax:',
x,
3181 & sqrt(-qsum2),
getpdf(
x,scalefacm*sqrt(-qsum2),
'GG'),
'gg',
3187 fmax=
log(-scalefacm**2*qsum2/q0**2)
3189 IF(qsum2.EQ.0.d0)
THEN
3196 write(logfid,*)
'x,sqrt(qsum^2),getpdf:,fmax',
x,
3197 & sqrt(-qsum2),
getpdf(
x,scalefacm*sqrt(-qsum2),
'GQ'),
'gq',
3203 write(logfid,*)
'x,sqrt(qsum^2),getpdf,fmax:',
x,
3204 & sqrt(-qsum2),
getpdf(
x,scalefacm*sqrt(-qsum2),
'QQ'),
'qq',
3212 &
write(logfid,*)
'error: weight=',
weight
3221 qsumvec(1)=allqs(nend,2)
3222 qsumvec(2)=allqs(nend,3)
3223 qsumvec(3)=allqs(nend,4)
3224 qsumvec(4)=allqs(nend,5)
3226 IF(-allqs(nend,1).GT.
low)
THEN
3231 deltat=allqs(nend,6)-tstart
3238 IF(((tleft.LT.tauest).OR.(
pyr(0).GT.1.d0/(nend*1.d0)))
3239 & .AND.(.NOT.norad))
THEN
3244 qsumvec(1)=allqs(nend,2)
3245 qsumvec(2)=allqs(nend,3)
3246 qsumvec(3)=allqs(nend,4)
3247 qsumvec(4)=allqs(nend,5)
3249 IF(-allqs(nend,1).GT.
low)
THEN
3254 deltat=allqs(nend,6)-tstart
3288 common/hepmcid/hpmcfid,logfid
3289 integer hpmcfid,logfid
3293 DOUBLE PRECISION p,
v
3295 common/time/mv(23000,5)
3298 common/ftimefac/ftfac
3299 DOUBLE PRECISION ftfac
3301 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
3304 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
3305 LOGICAL angord,scatrecoil,allhad,
compress
3307 common/disc/ndisc,nstrange,ngood,errcount,wdisc,discard
3309 INTEGER ndisc,nstrange,ngood,errcount
3310 double precision wdisc
3312 common/
coherent/nstart,nend,allqs(10000,6),scatcentres(10000,10),
3315 DOUBLE PRECISION allqs,scatcentres,qsumvec,qsum2
3317 INTEGER l,
counter,countmax,count2
3318 DOUBLE PRECISION tnow,deltat,newmass,tleft,deltal,q2old,
3326 IF(-qsum2.GT.
p(l,4)**2)
3327 &
write(logfid,*)
'DOFISTATESCAT has a problem:',-qsum2,
p(l,4)**2
3329 IF(
k(l,2).EQ.21)
THEN
3334 low=q0**2/scalefacm**2
3342 IF(overq0.OR.(-qsum2.GT.
low))
THEN
3343 newmass=
getnewmass(l,scalefacm**2*qsum2,scalefacm**2*q2old,
3344 & newmass,.
false.,1.d0,
z,qqbar)
3353 222
IF((newmass.EQ.0.d0).OR.(newmass.EQ.
p(l,5)))
THEN
3356 tauest=ftfac*(1.-
phi)*0.2*
p(l,4)/newmass**2
3359 tauest=-
log(1.d0-rtau)/lambda
3360 noscat=.NOT.
getdeltat(l,tnow+tsum,min(tauest,tleft),deltal)
3364 IF(nend.gt.countmax)
THEN
3368 IF(nstart.EQ.0) nstart=1
3371 IF((newmass.NE.0.d0).AND.(newmass.NE.
p(l,5)))
THEN
3372 phi=
phi+5.*deltal*newmass**2/(1.*
p(l,4))
3374 allqs(nend,6)=tnow+tsum
3375 qsumvecold(1)=qsumvec(1)
3376 qsumvecold(2)=qsumvec(2)
3377 qsumvecold(3)=qsumvec(3)
3378 qsumvecold(4)=qsumvec(4)
3382 219 CALL
getqvec(l,nend,tnow+tsum-mv(l,4),1.d0)
3383 IF(-qsum2.GT.
p(l,4)**2)
THEN
3384 qsumvec(1)=qsumvecold(1)
3385 qsumvec(2)=qsumvecold(2)
3386 qsumvec(3)=qsumvecold(3)
3387 qsumvec(4)=qsumvecold(4)
3389 IF(count2.LT.100)
THEN
3401 IF(overq0.OR.(-qsum2.GT.
low))
THEN
3402 newmass=
getnewmass(l,scalefacm**2*qsum2,scalefacm**2*q2old,
3403 & newmass,.
false.,1.d0,
z,qqbar)
3409 218
if ((newmass**2.gt.
low).and.(newmass.ne.
p(l,5)))
then
3410 if ((tleft.LT.tauest).OR.(
pyr(0).GT.1.d0/(nend*1.d0)))
then
3411 if (nend.eq.countmax)
then
3413 else if (tleft.LT.tauest)
then
3427 qsumvec(1)=allqs(nend,2)
3428 qsumvec(2)=allqs(nend,3)
3429 qsumvec(3)=allqs(nend,4)
3430 qsumvec(4)=allqs(nend,5)
3431 IF(-allqs(nend,1).GT.
low)
THEN
3451 DOUBLE PRECISION p,
v
3453 common/time/mv(23000,5)
3456 common/angor/za(23000),zd(23000),thetaa(23000),qqbard(23000)
3457 DOUBLE PRECISION za,zd,thetaa
3460 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
3463 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
3464 LOGICAL angord,scatrecoil,allhad,
compress
3467 DOUBLE PRECISION q2,qold2,
r,
pyr,pnosplit1,pnosplit2,
z,qa,
3469 LOGICAL in,qqbardec,qqbarold
3472 IF(
x*
p(l,4).LT.q0)
THEN
3478 IF (-q2.LT.q0**2)
THEN
3482 IF(
k(l,2).EQ.21)
THEN
3487 IF(sqrt(-qold2).LE.q0)
THEN
3490 &
x*
p(l,4),typ,
x*
p(l,4),
in,zdec,qqbardec)
3493 & sqrt(-q2),
in,zdec,qqbardec)
3501 IF(-q2.GT.-qold2)
THEN
3504 qtmp=
getmass(0.d0,sqrt(-q2),-1.d0,
x*
p(l,4),typ,
3505 & sqrt(-q2),
in,zdec,qqbardec)
3506 IF(qtmp.LT.sqrt(-qold2))
THEN
3518 pkeep=(1.-pnosplit2)/(1.-pnosplit1)
3519 IF(
pyr(0).LT.pkeep)
THEN
3520 IF(
p(l,5).LT.sqrt(-q2))
THEN
3524 & sqrt(-q2),
in,zdec,qqbardec)
3534 IF(-q2.GT.-qold2)
THEN
3536 &
x*
p(l,4),typ,
x*
p(l,4),
in,zdec,qqbardec)
3558 common/hepmcid/hpmcfid,logfid
3559 integer hpmcfid,logfid
3563 DOUBLE PRECISION p,
v
3565 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
3568 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
3569 LOGICAL angord,scatrecoil,allhad,
compress
3573 &
scatprimfunc,ms1,md1,shat,pcms2,avmom(5),
x,
y,
z,
t,
getmd
3577 &avmom(1),avmom(2),avmom(3),avmom(4),avmom(5))
3579 shat = avmom(5)**2 +
p(
line,5)**2 + 2.*(avmom(4)*
p(
line,4)
3581 pcms2 = (shat+
p(
line,5)**2-ms1**2)**2/(4.*shat)-
p(
line,5)**2
3583 low=q0**2/scalefacm**2
3584 IF((up.LE.
low).OR.(
p(
line,4).LT.q0/scalefacm))
THEN
3588 IF(
k(
line,2).EQ.21)
THEN
3593 IF(sigmatot.EQ.0.d0)
THEN
3605 IF(sigmatot.EQ.0.d0)
THEN
3626 common/hepmcid/hpmcfid,logfid
3627 integer hpmcfid,logfid
3631 DOUBLE PRECISION p,
v
3633 common/time/mv(23000,5)
3636 common/
coherent/nstart,nend,allqs(10000,6),scatcentres(10000,10),
3639 DOUBLE PRECISION allqs,scatcentres,qsumvec,qsum2
3641 common/disc/ndisc,nstrange,ngood,errcount,wdisc,discard
3643 INTEGER ndisc,nstrange,ngood,errcount
3644 double precision wdisc
3646 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
3649 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
3650 LOGICAL angord,scatrecoil,allhad,
compress
3653 DOUBLE PRECISION xsc,ysc,zsc,tsc,
getmd,
gettemp,
dt,
x,
pyr,newmom(4),
3654 &
t,
pt,maxt,phi2,
beta(3),
phi,theta,
gett,
pyp,
pi,
pt2,
getms,
3655 &savemom(5),theta2,mb2,
pz,
kt2,phiq,maxt2,xi,md,shat,pcms2,
3658 DATA pi/3.141592653589793d0/
3669 xsc=mv(l,1)+
dt*
p(l,1)/
p(l,4)
3670 ysc=mv(l,2)+
dt*
p(l,2)/
p(l,4)
3671 zsc=mv(l,3)+
dt*
p(l,3)/
p(l,4)
3673 md =
getmd(xsc,ysc,zsc,tsc)
3676 &avmom(1),avmom(2),avmom(3),avmom(4),avmom(5))
3682 xi = sqrt(
max(
x**2*
p(l,4)**2,
p(l,5)**2) -
p(l,5)**2)/
pyp(l,8)
3690 &
k(1,2),
p(1,1),
p(1,2),
p(1,3),
p(1,4),
p(1,5))
3696 IF(
k(1,2).EQ.21)typs=
'G'
3698 shat = avmom(5)**2 + savemom(5)**2 + 2.*(avmom(4)*savemom(4)
3699 & -avmom(1)*savemom(1)-avmom(2)*savemom(2)-avmom(3)*savemom(3))
3700 pcms2 = (shat+savemom(5)**2-avmom(5)**2)**2/(4.*shat)
3705 scatcentres(
j,1)=
k(1,2)
3706 scatcentres(
j,2)=
p(1,1)
3707 scatcentres(
j,3)=
p(1,2)
3708 scatcentres(
j,4)=
p(1,3)
3709 scatcentres(
j,5)=
p(1,4)
3710 scatcentres(
j,6)=
p(1,5)
3711 scatcentres(
j,7)=mv(1,1)
3712 scatcentres(
j,8)=mv(1,2)
3713 scatcentres(
j,9)=mv(1,3)
3714 scatcentres(
j,10)=mv(1,4)
3725 CALL
pyrobo(l,l,-theta,0d0,0d0,0d0,0d0)
3726 CALL
pyrobo(1,1,-theta,0d0,0d0,0d0,0d0)
3728 204
t=-
gett(0.d0,maxt,md)
3729 202 newmom(4)=
p(l,4)+
t/(2.*
p(1,5))
3730 newmom(3)=(
t-2.*
p(l,5)**2+2.*
p(l,4)*newmom(4))/(2.*
p(l,3))
3731 pt2=newmom(4)**2-newmom(3)**2-
p(l,5)**2
3732 IF(dabs(
pt2).LT.1.d-10)
pt2=0.d0
3733 IF(
t.EQ.0.d0)
pt2=0.d0
3740 newmom(1)=
pt*cos(phi2)
3741 newmom(2)=
pt*sin(phi2)
3742 p(1,1)=newmom(1)-
p(l,1)
3743 p(1,2)=newmom(2)-
p(l,2)
3744 p(1,3)=newmom(3)-
p(l,3)
3745 p(1,4)=newmom(4)-
p(l,4)
3748 CALL
pyrobo(l,l,theta,0d0,0d0,0d0,0d0)
3749 CALL
pyrobo(1,1,theta,0d0,0d0,0d0,0d0)
3759 qsumvec(1)=qsumvec(1)+allqs(nend,2)
3760 qsumvec(2)=qsumvec(2)+allqs(nend,3)
3761 qsumvec(3)=qsumvec(3)+allqs(nend,4)
3762 qsumvec(4)=qsumvec(4)+allqs(nend,5)
3763 qsum2=qsumvec(4)**2-qsumvec(1)**2-qsumvec(2)**2-qsumvec(3)**2
3764 IF(qsum2.GT.0.d0)
THEN
3765 qsumvec(1)=qsumvec(1)-allqs(nend,2)
3766 qsumvec(2)=qsumvec(2)-allqs(nend,3)
3767 qsumvec(3)=qsumvec(3)-allqs(nend,4)
3768 qsumvec(4)=qsumvec(4)-allqs(nend,5)
3769 qsum2=qsumvec(4)**2-qsumvec(1)**2-qsumvec(2)**2-qsumvec(3)**2
3771 write(logfid,*)
'GETQVEC unable to find q vector'
3794 common/hepmcid/hpmcfid,logfid
3795 integer hpmcfid,logfid
3799 DOUBLE PRECISION p,
v
3801 common/time/mv(23000,5)
3804 common/ftimefac/ftfac
3805 DOUBLE PRECISION ftfac
3807 common/colour/trip(23000),anti(23000),colmax
3808 INTEGER trip,anti,colmax
3810 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
3813 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
3814 LOGICAL angord,scatrecoil,allhad,
compress
3816 common/disc/ndisc,nstrange,ngood,errcount,wdisc,discard
3818 INTEGER ndisc,nstrange,ngood,errcount
3819 double precision wdisc
3821 common/angor/za(23000),zd(23000),thetaa(23000),qqbard(23000)
3822 DOUBLE PRECISION za,zd,thetaa
3825 common/
coherent/nstart,nend,allqs(10000,6),scatcentres(10000,10),
3828 DOUBLE PRECISION allqs,scatcentres,qsumvec,qsum2
3830 common/check/nscat,nscateff,nsplit
3831 DOUBLE PRECISION nscat,nscateff,nsplit
3833 common/
weight/evweight,sumofweights
3834 double precision evweight,sumofweights
3836 common/storescatcen/nscatcen,maxnscatcen,scatflav(10000),
3837 &scatcen(10000,5),writescatcen,writedummies
3838 integer nscatcen,maxnscatcen,scatflav
3839 double precision scatcen
3840 logical writescatcen,writedummies
3842 INTEGER l,
line,n1,n2,
j,
dir,lold,nold,colmaxold,statold,nscatcenold
3844 &newmass,deltam,dm,ttot,dmleft,lambda,time,endtime,
x,tmp,
3845 &m32,newm2,shat,theta2,
z,
gettemp,e3new,e4new,p32,p42,p3old,
3846 &newm,mass2,enew,
pt2,
pt,
pl,m12,firsttime,pcms2
3848 LOGICAL retrysplit,qqbar,qqbardec,rejectt,redokin,reshuffle
3849 DATA pi/3.141592653589793d0/
3851 IF((
n+2*(n2-n1+1)).GT.22990)
THEN
3852 write(logfid,*)
'event too long for event record'
3865 204 deltam=newm2-
p(l,5)
3870 ttot=ttot+allqs(
j,1)
3878 IF(
k(
line,2).EQ.21)
THEN
3880 IF(
pyr(0).LT.0.5)
THEN
3890 k(1,2)=scatcentres(
j,1)
3891 p(1,1)=scatcentres(
j,2)
3892 p(1,2)=scatcentres(
j,3)
3893 p(1,3)=scatcentres(
j,4)
3894 p(1,4)=scatcentres(
j,5)
3895 p(1,5)=scatcentres(
j,6)
3896 mv(1,1)=scatcentres(
j,7)
3897 mv(1,2)=scatcentres(
j,8)
3898 mv(1,3)=scatcentres(
j,9)
3899 mv(1,4)=scatcentres(
j,10)
3911 IF ((
beta(1).GT.1.d0).OR.(
beta(2).GT.1.d0).OR.(
beta(3).GT.1.d0)
3912 & .or.(sqrt(
beta(1)**2+
beta(2)**2+
beta(3)**2).gt.1.d0))
THEN
3917 205
if (.not.reshuffle)
then
3928 CALL
pyrobo(1,1,-theta,0d0,0d0,0d0,0d0)
3930 maxt = -2.*
p(1,5)*
p(
line,4)
3936 203 enew =
p(
line,4)+
t/(2.*
p(1,5))
3938 pt2 = enew**2-
pl**2-m12
3939 if (
t.eq.0.d0)
pt2 = 0.d0
3940 if (dabs(
pt2).lt.1.d-8)
pt2 = 0.d0
3941 if (
pt2.lt.0.d0)
then
3942 write(logfid,*)
' This should not have happened: pt^2<0!'
3943 write(logfid,*)
t,enew,
pl,
pt2
3961 mass2 =
p(
n-1,4)**2-
p(
n-1,1)**2-
p(
n-1,2)**2-
p(
n-1,3)**2
3962 if ((mass2.lt.0.d0).and.(mass2.gt.-1.-6)) mass2=0.d0
3964 &
write(logfid,*)
'messed up scattering centres mass^2: ',
3966 p(
n-1,5)=sqrt(mass2)
3967 if (abs(
p(
n-1,5)-
p(1,5)).gt.1.d-6)
3968 &
write(logfid,*)
'messed up scattering centres mass: ',
3969 &
p(
n-1,5),
p(1,5),
p(l,5)
3981 if ((
p(1,4).lt.0.d0).or.(
p(
line,4).lt.0.d0))
then
3992 CALL
pyrobo(1,1,-theta,0d0,0d0,0d0,0d0)
3993 shat = (
p(1,4)+
p(
line,4))**2
3996 maxt = -4.*
p(
line,3)**2
4001 theta2 = acos(1.d0+
t/(2.*
p(
line,3)**2))
4004 p(
n,1)=
p(
line,3)*sin(theta2)*cos(phi2)
4005 p(
n,2)=
p(
line,3)*sin(theta2)*sin(phi2)
4014 mass2 =
p(
n-1,4)**2-
p(
n-1,1)**2-
p(
n-1,2)**2-
p(
n-1,3)**2
4015 if ((mass2.lt.0.d0).and.(mass2.gt.-1.-6)) mass2=0.d0
4017 &
write(logfid,*)
'messed up scattering centres mass^2: ',
4019 p(
n-1,5)=sqrt(mass2)
4020 if (abs(
p(
n-1,5)-
p(1,5)).gt.1.d-6)
4021 &
write(logfid,*)
'messed up scattering centres mass: ',
4022 &
p(
n-1,5),
p(1,5),
p(l,5)
4042 IF(allhad.and.(.not.rejectt))
THEN
4043 IF(
k(
n,2).EQ.21)
THEN
4051 ELSEIF(
k(
n,2).GT.0)
THEN
4064 IF(
k(
line,1).EQ.1)
THEN
4082 IF(allhad.and.(.not.rejectt))
THEN
4083 IF((
k(
n,2).GT.0).AND.(
dir.GE.0))
THEN
4084 trip(
n-1)=trip(
line)
4088 anti(
n-1)=anti(
line)
4095 if (reshuffle.and.(dm.gt.0.d0))
then
4098 IF(ttot.EQ.0.d0)
THEN
4101 if (dmleft.lt.0.d0)
then
4102 dm=
max(dmleft*
t/ttot*1.5d0,dmleft)
4104 dm=min(dmleft*
t/ttot*1.5d0,dmleft)
4107 ttot=ttot-allqs(
j,1)
4110 if (newmass.lt.0.d0)
then
4115 e3new = (shat + m32 -
p(1,5)**2)/(2.d0*sqrt(shat))
4116 e4new = (shat - m32 +
p(1,5)**2)/(2.d0*sqrt(shat))
4117 p32 = e3new**2 - m32
4118 p42 = e4new**2 -
p(1,5)**2
4119 if ((p32.lt.0.d0).or.(p42.lt.0.d0).or.
4120 & (e3new.lt.0.d0).or.(e4new.lt.0.d0))
then
4124 e3new = sqrt(shat) - e4new
4126 if ((e3new.lt.0.d0).or.(e4new.lt.0.d0))
then
4131 if (
p(
n,5).lt.0.d0)
then
4138 p(
n,1) = sqrt(p32)*
p(
n,1)/p3old
4139 p(
n,2) = sqrt(p32)*
p(
n,2)/p3old
4140 p(
n,3) = sqrt(p32)*
p(
n,3)/p3old
4142 p(
n,5) =
sign(sqrt(abs(m32)),newmass)
4143 tmp =
p(
n,4)**2-
p(
n,1)**2-
p(
n,2)**2-
p(
n,3)**2
4144 if (abs(tmp-m32).gt.1.d-6)
4145 &
write(logfid,*)
'Oups, messed up projectiles mass:',
4148 p(
n-1,1) = sqrt(p42)*
p(
n-1,1)/p3old
4149 p(
n-1,2) = sqrt(p42)*
p(
n-1,2)/p3old
4150 p(
n-1,3) = sqrt(p42)*
p(
n-1,3)/p3old
4152 tmp =
p(
n-1,4)**2-
p(
n-1,1)**2-
p(
n-1,2)**2-
p(
n-1,3)**2
4154 if (abs(tmp).gt.1.d-6)
4155 &
write(logfid,*)
'Oups, messed up scattering centres mass:',
4156 & tmp,p3old,
p(
n-1,1),
p(
n-1,2),
p(
n-1,3),
p(
n-1,4),
p(
n-1,5)
4157 if ((abs(
p(
n,1)+
p(
n-1,1)).gt.1.d-6).or.
4158 & (abs(
p(
n,2)+
p(
n-1,2)).gt.1.d-6).or.
4159 & (abs(
p(
n,3)+
p(
n-1,3)).gt.1.d-6))
4160 &
write(logfid,*)
'Oups, momentum not conserved',
4172 CALL
pyrobo(
n-1,
n,theta,0d0,0d0,0d0,0d0)
4178 CALL
pyrobo(1,1,theta,0d0,0d0,0d0,0d0)
4181 if (.not.allhad)
then
4184 IF(scatrecoil.AND.(
p(
n-1,4).GT.(10.*3.*
4185 &
gettemp(mv(1,1),mv(1,2),mv(1,3),mv(1,4)))))
THEN
4191 if (rejectt)
k(
n-1,1)=11
4195 mv(
n-1,1)=mv(
line,1)
4197 mv(
n-1,2)=mv(
line,2)
4199 mv(
n-1,3)=mv(
line,3)
4207 IF(
p(
n-1,5).GT.
p(1,5))
THEN
4208 lambda=1.d0/(ftfac*0.2*
p(
n-1,4)/
p(
n-1,5)**2)
4209 mv(
n-1,5)=mv(
n-1,4)-
log(1.d0-
pyr(0))/lambda
4214 mv(
n,5)=scatcentres(
j+1,10)
4216 IF(
p(
n,5).GT.0.d0)
THEN
4217 IF(deltam.EQ.0.d0)
THEN
4221 lambda=1.d0/(ftfac*
p(
n,4)*0.2/
p(
n,5)**2)
4222 endtime=scatcentres(
j,10)-
log(1.d0-
pyr(0))/lambda
4232 mv(
line,5)=allqs(
j,6)
4236 if (writescatcen.and.(.not.rejectt).and.
4237 & (nscatcen.lt.maxnscatcen))
then
4238 nscatcen = nscatcen+1
4239 if (nscatcen.le.maxnscatcen)
then
4240 scatflav(nscatcen) =
k(1,2)
4241 scatcen(nscatcen,1) =
p(1,1)
4242 scatcen(nscatcen,2) =
p(1,2)
4243 scatcen(nscatcen,3) =
p(1,3)
4244 scatcen(nscatcen,4) =
p(1,4)
4245 scatcen(nscatcen,5) =
p(1,5)
4248 &
'WARNING: no room left to store further scattering centres'
4262 dmleft=dmleft-(
p(
n,5)-
p(
line,5))
4264 tmp = abs(
p(
n,4)**2-
p(
n,1)**2-
p(
n,2)**2-
p(
n,3)**2)-
p(
n,5)**2
4265 if (abs(tmp).ge.1.d-6)
4266 &
write(logfid,*)tmp,
j,
p(l,5),
p(
line,5),
p(
n,5)
4268 if (
p(
n,5).lt.0.d0)
then
4272 if (
p(
n,5).ne.newm2)
then
4278 if (
p(l,5).le.0.d0)
then
4281 if (
p(l,5).lt.q0)
then
4282 if ((newm2.eq.newm).and.(newm.ne.q0+1.d-6))
then
4295 if ((
k(
n,1).eq.1).and.
4296 & ((
p(
n,5).lt.0.d0).or.((
p(
n,5).gt.0.d0).and.(
p(
n,5).lt.q0))))
4297 &
write(logfid,*)
'dokinematics did not reach sensible mass: ',
4298 &
p(
n,5),newm,
p(l,5),newm2
4299 nscateff=nscateff+evweight
4307 DOUBLE PRECISION FUNCTION getproba(QI,QF,QAA,ZAA,EBB,TYPE,
4311 common/sudaint/qa,za2,eb,
t,instate,typ
4312 DOUBLE PRECISION qa,za2,eb,
t
4338 common/hepmcid/hpmcfid,logfid
4339 integer hpmcfid,logfid
4341 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
4344 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
4345 LOGICAL angord,scatrecoil,allhad,
compress
4347 common/sudaint/qa,za2,eb,
t,instate,typ
4348 DOUBLE PRECISION qa,za2,eb,
t
4352 DOUBLE PRECISION qmax1,qa1,qb1,za1,eb1,tmax,tb,ystart,epsi,
4360 IF(qb2.LT.q0)
write(logfid,*)
'error: Q < Q0',qb2,qmax1
4361 IF(qb2.LT.(q0+1.d-10)) qb2=qb2+1.d-10
4363 IF(qb2.LT.q0)
write(logfid,*)
'error: Q < min',qb2,qmax1
4364 IF(qb2.LT.(q0+1.d-10)) qb2=qb2+1.d-10
4366 IF(qb2.GE.(qmax1-1.d-10))
THEN
4378 hfirst=0.01*(qmax1-qb1)
4380 CALL
odeint(ystart,qb2,qmax1,epsi,hfirst,0.d0,1)
4393 common/hepmcid/hpmcfid,logfid
4394 integer hpmcfid,logfid
4396 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
4399 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
4400 LOGICAL angord,scatrecoil,allhad,
compress
4402 common/sudaint/qa,za2,eb,
t,instate,typ
4403 DOUBLE PRECISION qa,za2,eb,
t
4407 DOUBLE PRECISION qmax1,qb,qb1,za1,ea1,ystart,epsi,
4413 IF(qb1.LT.q0)
write(logfid,*)
'error: Q < Q0',qb1,qmax1
4414 IF(qb1.LT.(q0+1.d-12)) qb1=qb1+1.d-12
4415 IF(qb1.GE.(qmax1-1.d-12))
THEN
4419 hfirst=0.01*(qmax1-qb1)
4421 CALL
odeint(ystart,qb1,qmax1,epsi,hfirst,0.d0,6)
4430 DOUBLE PRECISION FUNCTION deriv(XVAL,W4)
4433 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
4436 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
4437 LOGICAL angord,scatrecoil,allhad,
compress
4439 common/intsplitf/qquad,fm
4440 DOUBLE PRECISION qquad,fm
4442 common/sudaint/qa,za2,eb,
t,instate,typ
4443 DOUBLE PRECISION qa,za2,eb,
t
4447 common/pdfintv/
xmax,
z
4448 DOUBLE PRECISION xmax,
z
4450 common/xsecv/qlow,mdx
4451 DOUBLE PRECISION qlow,mdx
4457 DATA pi/3.141592653589793d0/
4468 deriv=(1.+fm)*
alphas(xval*(1.-xval)*qquad/1.,lps)*
4473 & *
pgg(xval)/(2.*
pi)
4479 deriv=exp(-xval)/xval
4484 & *
alphas((1.-
z)*xval**2/1.,lps)
4488 & *
alphas((1.-
z)*xval**2/1.,lps)
4492 & *
alphas((1.-
z)*xval**2/1.,lps)
4494 ELSEIF(w4.EQ.10)
THEN
4496 & *
alphas((1.-
z)*xval**2/1.,lps)*
4497 & *2.*
pgg(
z)/(2.*
pi*xval)
4498 ELSEIF(w4.EQ.11)
THEN
4501 ELSEIF(w4.EQ.12)
THEN
4504 ELSEIF(w4.EQ.13)
THEN
4506 & *3.*2.*
pi*
alphas(xval+mdx**2,lqcd)**2/(2.*(xval+mdx**2)**2)
4507 ELSEIF(w4.EQ.14)
THEN
4509 & *2.*2.*
pi*
alphas(xval+mdx**2,lqcd)**2/(3.*(xval+mdx**2)**2)
4510 ELSEIF(w4.EQ.21)
THEN
4513 ELSEIF(w4.EQ.22)
THEN
4516 ELSEIF(w4.EQ.23)
THEN
4519 ELSEIF(w4.EQ.24)
THEN
4534 common/hepmcid/hpmcfid,logfid
4535 integer hpmcfid,logfid
4537 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
4540 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
4541 LOGICAL angord,scatrecoil,allhad,
compress
4543 common/splitint/splitiggv(1000,1000),splitiqqv(1000,1000),
4544 &splitiqgv(1000,1000),qval(1000),zmval(1000),qmax,zmmin,npoint
4546 DOUBLE PRECISION splitiggv,splitiqqv,splitiqgv,
4547 &qval,zmval,qmax,zmmin
4549 common/intsplitf/qquad,fm
4550 DOUBLE PRECISION qquad,fm
4552 common/extrapolations/ntotspliti,noverspliti,ntotpdf,noverpdf,
4553 &ntotxsec,noverxsec,ntotsuda,noversuda
4554 integer ntotspliti,noverspliti,ntotpdf,noverpdf,
4555 &ntotxsec,noverxsec,ntotsuda,noversuda
4557 INTEGER i,
j,lt,qlmax,zlmax,qline,zline
4558 DOUBLE PRECISION qa,qb,zeta,eb,
low,x1a(2),x2a(2),ya(2,2),
y,
4559 &splitintgg,splitintqg,a,b,yb(2)
4562 ntotspliti=ntotspliti+1
4563 if (qb.gt.qmax)
then
4564 noverspliti=noverspliti+1
4565 if (noverspliti.le.25)
4566 &
write(logfid,*)
'WARNING in getspliti: need to extrapolate: ',
4571 IF(angord.AND.(zeta.NE.1.d0))
THEN
4572 low=
max(0.5-0.5*sqrt(1.-q0**2/qb**2)
4573 & *sqrt(1.-qb**2/eb**2),
4574 & 0.5-0.5*sqrt(1.-4.*qb**2*(1.-zeta)/(zeta*qa**2)))
4576 low=0.5-0.5*sqrt(1.-q0**2/qb**2)
4577 & *sqrt(1.-qb**2/eb**2)
4580 qlmax=int((qb-qval(1))*npoint/(qval(1000)-qval(1))+1)
4582 qline=min(qline,npoint)
4583 zlmax=int((
log(
low)-
log(zmval(1)))*npoint/
4584 & (
log(zmval(1000))-
log(zmval(1)))+1)
4586 zline=min(zline,npoint)
4587 IF((qline.GT.999).OR.(zline.GT.999).OR.
4588 & (qline.LT.1).OR.(zline.LT.1))
THEN
4589 write(logfid,*)
'ERROR in GETSPLITI: line number out of bound',
4592 IF((type1.EQ.
'GG').OR.(type1.EQ.
'GC'))
THEN
4594 x1a(
i)=qval(qline-1+
i)
4595 x2a(
i)=zmval(zline-1+
i)
4597 ya(
i,
j)=splitiggv(qline-1+
i,zline-1+
j)
4601 a=(ya(
i,2)-ya(
i,1))/(x2a(2)-x2a(1))
4605 IF(x1a(1).EQ.x1a(2))
THEN
4608 a=(yb(2)-yb(1))/(x1a(2)-x1a(1))
4612 IF(type1.EQ.
'GG')
THEN
4615 splitintgg=min(
y,10.d0)
4618 IF((type1.EQ.
'QG').OR.(type1.EQ.
'GC'))
THEN
4620 x1a(
i)=qval(qline-1+
i)
4621 x2a(
i)=zmval(zline-1+
i)
4623 ya(
i,
j)=splitiqgv(qline-1+
i,zline-1+
j)
4627 a=(ya(
i,2)-ya(
i,1))/(x2a(2)-x2a(1))
4631 IF(x1a(1).EQ.x1a(2))
THEN
4634 a=(yb(2)-yb(1))/(x1a(2)-x1a(1))
4638 IF(type1.EQ.
'QG')
THEN
4641 splitintqg=nf*min(
y,10.d0)
4644 IF(type1.EQ.
'QQ')
THEN
4646 x1a(
i)=qval(qline-1+
i)
4647 x2a(
i)=zmval(zline-1+
i)
4649 ya(
i,
j)=splitiqqv(qline-1+
i,zline-1+
j)
4653 a=(ya(
i,2)-ya(
i,1))/(x2a(2)-x2a(1))
4657 IF(x1a(1).EQ.x1a(2))
THEN
4660 a=(yb(2)-yb(1))/(x1a(2)-x1a(1))
4666 IF(type1.EQ.
'GC')
getspliti=splitintgg+splitintqg
4676 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
4679 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
4680 LOGICAL angord,scatrecoil,allhad,
compress
4682 DOUBLE PRECISION qb,
low,
pi,
y,splitintgg,splitintqg,up,
ei
4684 DATA pi/3.141592653589793d0/
4687 up = 1. - q0**2/(4.*qb**2)
4688 IF((type1.EQ.
'GG').OR.(type1.EQ.
'GC'))
THEN
4695 & - lps**2*
ei(
log((1.-
low)*qb**2/lps**2))/qb**2
4696 & + lps**4*
ei(2.*
log((1.-
low)*qb**2/lps**2))/qb**4
4697 & - lps**6*
ei(3.*
log((1.-
low)*qb**2/lps**2))/qb**6
4698 & -
log(
log((1.-up)*qb**2/lps**2))
4699 & + lps**2*
ei(
log((1.-up)*qb**2/lps**2))/qb**2
4700 & - lps**4*
ei(2.*
log((1.-up)*qb**2/lps**2))/qb**4
4701 & + lps**6*
ei(3.*
log((1.-up)*qb**2/lps**2))/qb**6
4703 & *3.*12.*
pi/(2.*
pi*(33.-2.*nf))
4704 IF(type1.EQ.
'GG')
THEN
4710 IF((type1.EQ.
'QG').OR.(type1.EQ.
'GC'))
THEN
4716 y = ( 2.*lps**6*
ei(3.*
log((1.-
low)*qb**2/lps**2))/qb**6
4717 & - 2.*lps**4*
ei(2.*
log((1.-
low)*qb**2/lps**2))/qb**4
4718 & + 2.*lps**2*
ei(
log((1.-
low)*qb**2/lps**2))/qb**2
4719 & - 2.*lps**6*
ei(3.*
log((1.-up)*qb**2/lps**2))/qb**6
4720 & + 2.*lps**4*
ei(2.*
log((1.-up)*qb**2/lps**2))/qb**4
4721 & - 2.*lps**2*
ei(
log((1.-up)*qb**2/lps**2))/qb**2 )
4722 & *12.*
pi/(2.*2.*
pi*(33.-2.*nf))
4723 IF(type1.EQ.
'QG')
THEN
4729 IF(type1.EQ.
'QQ')
THEN
4736 & - 2.*lps**2*
ei(
log((1.-
low)*qb**2/lps**2))/qb**2
4737 & + lps**4*
ei(2.*
log((1.-
low)*qb**2/lps**2))/qb**4
4738 & - 2.*
log(
log((1.-up)*qb**2/lps**2))
4739 & + 2.*lps**2*
ei(
log((1.-up)*qb**2/lps**2))/qb**2
4740 & - lps**4*
ei(2.*
log((1.-up)*qb**2/lps**2))/qb**4 )
4741 & *4.*12.*
pi/(3.*2.*
pi*(33.-2.*nf))
4744 IF(type1.EQ.
'GQ')
THEN
4751 & *4.*12.*
pi/(3.*2.*
pi*(33.-2.*nf)*
log(qb**2/lps**2))
4754 IF(type1.EQ.
'GC')
getinspliti=splitintgg+splitintqg
4764 common/hepmcid/hpmcfid,logfid
4765 integer hpmcfid,logfid
4767 common/pdfs/qinqx(2,1000),ginqx(2,1000),qingx(2,1000),
4769 DOUBLE PRECISION qinqx,ginqx,qingx,gingx
4771 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
4774 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
4775 LOGICAL angord,scatrecoil,allhad,
compress
4777 common/pdfintv/
xmax,
z
4778 DOUBLE PRECISION xmax,
z
4780 DOUBLE PRECISION x,
q,qlow,qhigh,ystart,epsi,hfirst
4784 IF((
x.LT.0.d0).OR.(
x.GT.1.d0).OR.(
q.LT.q0))
THEN
4785 write(logfid,*)
'error in GETPDF: parameter out of bound',
x,
q
4794 qlow=
max(q0,q0/(2.*sqrt(1.-
x)))
4796 IF((qlow.GE.qhigh*(1.d0-1.d-10)).OR.(
x.GT.1.d0-1.d-10))
THEN
4799 hfirst=0.01*(qhigh-qlow)
4801 CALL
odeint(ystart,qlow,qhigh,epsi,hfirst,0.d0,7)
4804 ELSEIF(typ.EQ.
'GQ')
THEN
4808 qlow=
max(q0,
max(q0/(2.*sqrt(
x)),q0/(2.*sqrt(1.-
x))))
4810 IF((qlow.GE.qhigh*(1.d0-1.d-10)).OR.(
x.LT.0.d0+1.d-10)
4811 & .OR.(
x.GT.1.d0-1.d-10))
THEN
4814 hfirst=0.01*(qhigh-qlow)
4816 CALL
odeint(ystart,qlow,qhigh,epsi,hfirst,0.d0,8)
4819 ELSEIF(typ.EQ.
'QG')
THEN
4823 qlow=
max(q0,q0/(2.*sqrt(1.-
x)))
4825 IF((qlow.GE.qhigh*(1.d0-1.d-10)).OR.(
x.GT.1.d0-1.d-10))
THEN
4828 hfirst=0.01*(qhigh-qlow)
4830 CALL
odeint(ystart,qlow,qhigh,epsi,hfirst,0.d0,9)
4833 ELSEIF(typ.EQ.
'GG')
THEN
4837 qlow=
max(q0,
max(q0/(2.*sqrt(
x)),q0/(2.*sqrt(1.-
x))))
4839 IF((qlow.GE.qhigh*(1.d0-1.d-10)).OR.(
x.LT.0.d0+1.d-10)
4840 & .OR.(
x.GT.1.d0-1d-10))
THEN
4843 hfirst=0.01*(qhigh-qlow)
4845 CALL
odeint(ystart,qlow,qhigh,epsi,hfirst,0.d0,10)
4849 write(logfid,*)
'error: pdf-type ',typ,
' does not exist'
4860 common/hepmcid/hpmcfid,logfid
4861 integer hpmcfid,logfid
4863 common/pdfs/qinqx(2,1000),ginqx(2,1000),qingx(2,1000),
4865 DOUBLE PRECISION qinqx,ginqx,qingx,gingx
4867 common/extrapolations/ntotspliti,noverspliti,ntotpdf,noverpdf,
4868 &ntotxsec,noverxsec,ntotsuda,noversuda
4869 integer ntotspliti,noverspliti,ntotpdf,noverpdf,
4870 &ntotxsec,noverxsec,ntotsuda,noversuda
4872 INTEGER j,q2close,q2line
4873 DOUBLE PRECISION q,xa(2),ya(2),
y,a,b
4877 if (
q**2.gt.qinqx(1,1000))
then
4880 &
write(logfid,*)
'WARNING in getpdfxint: need to extrapolate: ',
4881 &
q**2,qinqx(1,1000)
4884 q2close=int((
log(
q**2)-
log(qinqx(1,1)))*999.d0/
4885 & (
log(qinqx(1,1000))-
log(qinqx(1,1)))+1)
4886 q2line=
max(q2close,1)
4887 q2line=min(q2line,999)
4888 IF((q2line.GT.999).OR.(q2line.LT.1))
THEN
4889 write(logfid,*)
'ERROR in GETPDFXINT: line number out of bound',
4895 xa(
j)=qinqx(1,q2line-1+
j)
4896 ya(
j)=qinqx(2,q2line-1+
j)
4898 ELSEIF(typ.EQ.
'GQ')
THEN
4900 xa(
j)=ginqx(1,q2line-1+
j)
4901 ya(
j)=ginqx(2,q2line-1+
j)
4903 ELSEIF(typ.EQ.
'QG')
THEN
4905 xa(
j)=qingx(1,q2line-1+
j)
4906 ya(
j)=qingx(2,q2line-1+
j)
4908 ELSEIF(typ.EQ.
'GG')
THEN
4910 xa(
j)=gingx(1,q2line-1+
j)
4911 ya(
j)=gingx(2,q2line-1+
j)
4914 write(logfid,*)
'error in GETPDFXINT: unknown integral type ',typ
4916 a=(ya(2)-ya(1))/(xa(2)-xa(1))
4929 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
4932 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
4933 LOGICAL angord,scatrecoil,allhad,
compress
4935 common/pdfintv/
xmax,
z
4936 DOUBLE PRECISION xmax,
z
4938 DOUBLE PRECISION q,epsi,ystart,hfirst
4947 CALL
odeint(ystart,q0,
q,epsi,hfirst,0.d0,21)
4948 ELSEIF(typ.EQ.
'QG')
THEN
4949 CALL
odeint(ystart,q0,
q,epsi,hfirst,0.d0,23)
4950 ELSEIF(typ.EQ.
'GQ')
THEN
4951 CALL
odeint(ystart,q0,
q,epsi,hfirst,0.d0,22)
4952 ELSEIF(typ.EQ.
'GG')
THEN
4953 CALL
odeint(ystart,q0,
q,epsi,hfirst,0.d0,24)
4965 common/hepmcid/hpmcfid,logfid
4966 integer hpmcfid,logfid
4968 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
4971 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
4972 LOGICAL angord,scatrecoil,allhad,
compress
4974 common/xsecs/intq1(1001,101),intq2(1001,101),
4975 &intg1(1001,101),intg2(1001,101)
4976 DOUBLE PRECISION intq1,intq2,intg1,intg2
4978 common/xsecv/qlow,mdx
4979 DOUBLE PRECISION qlow,mdx
4981 common/extrapolations/ntotspliti,noverspliti,ntotpdf,noverpdf,
4982 &ntotxsec,noverxsec,ntotsuda,noversuda
4983 integer ntotspliti,noverspliti,ntotpdf,noverpdf,
4984 &ntotxsec,noverxsec,ntotsuda,noversuda
4986 INTEGER tline,tclose,mdclose,mdline,
i,
j
4987 DOUBLE PRECISION tm,x1a(2),x2a(2),ya(2,2),
y,md,yb(2),a,b
4991 if (tm.gt.intq1(1000,101))
then
4992 noverxsec=noverxsec+1
4994 &
write(logfid,*)
'WARNING in getxsecint: need to extrapolate: ',
4995 & tm,intq1(1000,101)
4998 tclose=int((
log(tm)-
log(intq1(1,101)))*999.d0/
4999 & (
log(intq1(1000,101))-
log(intq1(1,101)))+1)
5002 mdclose=int((md-intq1(1001,1))*99.d0/
5003 &(intq1(1001,100)-intq1(1001,1))+1)
5004 mdline=
max(mdclose,1)
5005 mdline=min(mdline,99)
5006 IF((
tline.GT.999).OR.(mdline.GT.99)
5007 & .OR.(
tline.LT.1).OR.(mdline.LT.1))
THEN
5008 write(logfid,*)
'ERROR in GETXSECINT: line number out of bound',
5012 IF(typ2.EQ.
'QA')
THEN
5015 x1a(
i)=intq1(1001,mdline-1+
i)
5021 ELSEIF(typ2.EQ.
'QB')
THEN
5024 x1a(
i)=intq2(1001,mdline-1+
i)
5030 ELSEIF(typ2.EQ.
'GA')
THEN
5033 x1a(
i)=intg1(1001,mdline-1+
i)
5039 ELSEIF(typ2.EQ.
'GB')
THEN
5042 x1a(
i)=intg2(1001,mdline-1+
i)
5049 write(logfid,*)
'error in GETXSECINT: unknown integral type ',
5053 a=(ya(
i,2)-ya(
i,1))/(x2a(2)-x2a(1))
5057 IF(x1a(1).EQ.x1a(2))
THEN
5060 a=(yb(2)-yb(1))/(x1a(2)-x1a(1))
5074 common/hepmcid/hpmcfid,logfid
5075 integer hpmcfid,logfid
5077 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
5080 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
5081 LOGICAL angord,scatrecoil,allhad,
compress
5088 ELSEIF(q1.LE.q0)
THEN
5095 write(logfid,*)
'ERROR: GETINSUDAFAST < 0:',
5108 common/hepmcid/hpmcfid,logfid
5109 integer hpmcfid,logfid
5111 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
5114 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
5115 LOGICAL angord,scatrecoil,allhad,
compress
5117 common/insuda/sudaqq(1000,2),sudaqg(1000,2),sudagg(1000,2),
5119 DOUBLE PRECISION sudaqq,sudaqg,sudagg,sudagc
5121 common/extrapolations/ntotspliti,noverspliti,ntotpdf,noverpdf,
5122 &ntotxsec,noverxsec,ntotsuda,noversuda
5123 integer ntotspliti,noverspliti,ntotpdf,noverpdf,
5124 &ntotxsec,noverxsec,ntotsuda,noversuda
5126 INTEGER qclose,qbin,
i
5127 DOUBLE PRECISION q,xa(2),ya(2),
y,a,b
5131 if (
q.gt.sudaqq(1000,1))
then
5132 noversuda=noversuda+1
5133 if (noversuda.le.25)
5134 &
write(logfid,*)
'WARNING in getinsudared: need to extrapolate: ',
5138 qclose=int((
log(
q)-
log(sudaqq(1,1)))*999.d0
5139 & /(
log(sudaqq(1000,1))-
log(sudaqq(1,1)))+1)
5142 IF((qbin.GT.999).OR.(qbin.LT.1))
THEN
5144 &
'ERROR in GETINSUDARED: line number out of bound',qbin
5146 IF(typ2.EQ.
'QQ')
THEN
5148 xa(
i)=sudaqq(qbin-1+
i,1)
5149 ya(
i)=sudaqq(qbin-1+
i,2)
5151 ELSEIF(typ2.EQ.
'QG')
THEN
5153 xa(
i)=sudaqg(qbin-1+
i,1)
5154 ya(
i)=sudaqg(qbin-1+
i,2)
5156 ELSEIF(typ2.EQ.
'GG')
THEN
5158 xa(
i)=sudagg(qbin-1+
i,1)
5159 ya(
i)=sudagg(qbin-1+
i,2)
5161 ELSEIF(typ2.EQ.
'GC')
THEN
5163 xa(
i)=sudagc(qbin-1+
i,1)
5164 ya(
i)=sudagc(qbin-1+
i,2)
5167 write(logfid,*)
'error in GETINSUDARED: unknown type ',typ2
5169 a=(ya(2)-ya(1))/(xa(2)-xa(1))
5174 write(logfid,*)
'ERROR: GETINSUDARED < 0:',
getinsudared,
q,typ2
5183 DOUBLE PRECISION FUNCTION getsscat(EN,px,py,PZ,MP,LW,TYPE1,TYPE2,
5187 common/hepmcid/hpmcfid,logfid
5188 integer hpmcfid,logfid
5190 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
5193 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
5194 LOGICAL angord,scatrecoil,allhad,
compress
5196 common/xsecv/qlow,mdx
5197 DOUBLE PRECISION qlow,mdx
5202 &
x,
y,
z,
t,
getmd,avmom(5),
px,
py,
getmdmin,
getmdmax,pproj,psct
5203 CHARACTER type1,type2
5205 IF(type1.EQ.
'Q')
THEN
5213 & avmom(1),avmom(2),avmom(3),avmom(4),avmom(5))
5214 shat = avmom(5)**2 + mp**2 +
5215 & 2.*(avmom(4)*en - avmom(1)*
px - avmom(2)*
py - avmom(3)*
pz)
5216 pcms2 = (shat+mp**2-avmom(5)**2)**2/(4.*shat)-mp**2
5224 call
maxscatcen(avmom(1),avmom(2),avmom(3),avmom(4),avmom(5))
5225 psct = sqrt(avmom(1)**2+avmom(2)**2+avmom(3)**2)
5226 pproj = sqrt(
px**2+
py**2+
pz**2)
5227 shat = avmom(5)**2 + mp**2 + 2.*(en*avmom(4) + pproj*psct)
5228 pcms2 = (shat+mp**2-avmom(5)**2)**2/(4.*shat)-mp**2
5236 IF((type2.EQ.
'C').OR.
5237 & ((type1.EQ.
'Q').AND.(type2.EQ.
'Q')).OR.
5238 & ((type1.EQ.
'G').AND.(type2.EQ.
'G')))
THEN
5243 low=q0**2/scalefacm**2
5245 IF(type1.EQ.
'Q')
THEN
5246 IF((type2.EQ.
'C').OR.(type2.EQ.
'G'))
THEN
5252 IF((type2.EQ.
'C').OR.(type2.EQ.
'G'))
THEN
5257 IF((type2.EQ.
'C').OR.(type2.EQ.
'Q'))
THEN
5265 &
write(logfid,*)
'error: cross section < 0',
getsscat,
'for',
5266 & en,mp,lw,type1,type2,lw**2,up
5275 DOUBLE PRECISION FUNCTION getmass(QBMIN,QBMAX,THETA,EP,TYPE,
5276 & max2,ins,zdec,qqbardec)
5279 common/hepmcid/hpmcfid,logfid
5280 integer hpmcfid,logfid
5284 DOUBLE PRECISION p,
v
5285 common/
pydat1/mstu(200),paru(200),mstj(200),parj(200)
5287 DOUBLE PRECISION paru,parj
5288 common/
pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
5289 INTEGER mdcy,mdme,kfdp
5290 DOUBLE PRECISION brat
5292 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
5295 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
5296 LOGICAL angord,scatrecoil,allhad,
compress
5298 common/time/mv(23000,5)
5301 common/alphasfac/ptfac
5302 DOUBLE PRECISION ptfac
5304 DOUBLE PRECISION qbmin,qbmax,theta,ep,max2,zdec,
5305 &q2min,alphmax,
alphas,log14,pref,q2max,sudaover,
gmin,
5307 &
r,
pyr,
z,rz,thetanew,
r2,
pi,
pqq,
pgg,
pqg,rmin
5309 LOGICAL ins,qqbardec
5310 DATA pi/3.141592653589793d0/
5314 alphmax =
alphas(3.*ptfac*q2min/16.,lps)
5317 IF(type.EQ.
'QQ')
THEN
5318 pref=4.*alphmax/(3.*2.*
pi)
5320 pref=29.*alphmax/(8.*2.*
pi)
5324 IF((qbmax.LE.qbmin).OR.(ep.LT.qbmin))
THEN
5334 21
if (q2max-qbmin**2.lt.1
e-4)
then
5337 IF(type.EQ.
'QQ')
THEN
5348 gmax = pref*
log(q2min/(4.*q2max))**2
5349 if (qbmin.gt.0.d0)
then
5350 rmin = exp(pref*
log(q2min/(4.*qbmin**2))**2-
gmax)
5355 r=
pyr(0)*(1.d0-rmin)+rmin
5367 cand = q2min*exp(sqrt(arg/pref))/4.
5368 eps = q2min/(4.*cand)
5370 if ((cand.lt.q2min).or.(cand.lt.qbmin**2))
then
5377 IF((cand.GT.max2**2).OR.(cand.GT.ep**2))
THEN
5388 trueeps=0.5-0.5*sqrt(1.-q2min/cand)
5389 & *sqrt(1.-cand/ep**2)
5391 &
WRITE(logfid,*)
'error in getmass: true eps < eps',trueeps,
eps
5394 if ((
z.lt.trueeps).or.(
z.gt.(1.-trueeps)))
then
5397 if (type.eq.
'QQ')
then
5403 oest = 2.*pref/(1.-
z)
5406 if (
pyr(0).lt.(17./29.))
z = 1.-
z
5411 trueval =
alphas(ptfac*
z*(1.-
z)*cand,lps)
5414 oest = alphmax*(17./(4.*
z)+3./(1.-
z))/(2.*
pi)
5417 thetanew = sqrt(cand/(
z*(1.-
z)))/ep
5418 if (angord.and.(theta.gt.0.).and.(thetanew.gt.theta))
5422 IF (
weight.GT.1.d0)
WRITE(logfid,*)
5423 &
'problem in getmass: weight> 1',
5433 IF(type.EQ.
'QQ')
THEN
5454 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
5457 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
5458 LOGICAL angord,scatrecoil,allhad,
compress
5460 DOUBLE PRECISION ti,ea,
eps,
pyr,
x,
r,help,
r1,epsi
5466 eps=
max(0.5-0.5*sqrt(1.-q0**2/ti)
5467 & *sqrt(1.-ti/ea**2),epsi)
5474 IF(type.EQ.
'QQ')
THEN
5477 IF(
r.LT.((1.+
x**2)/2.))
THEN
5482 ELSEIF(type.EQ.
'GG')
THEN
5483 x=1./(1.+((1.-
eps)/
eps)**(1.-2.*
r))
5485 help=((1.-
x)/
x+
x/(1.-
x)+
x*(1.-
x))/(1./(1.-
x)+1./
x)
5494 help=0.5*(
r**2+(1.-
r)**2)
5511 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
5514 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
5515 LOGICAL angord,scatrecoil,allhad,
compress
5518 DATA pi/3.141592653589793d0/
5521 & -
ei(-
log((
t+mdeb**2)/lqcd**2))/lqcd**2
5522 & - 1./((
t+mdeb**2)*
log((
t+mdeb**2)/lqcd**2)))/(33.-2.*nf)**2
5533 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
5536 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
5537 LOGICAL angord,scatrecoil,allhad,
compress
5539 DOUBLE PRECISION z,
q
5542 & +
log(1.-
z)))/((33.-2.*nf)*3.)
5553 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
5556 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
5557 LOGICAL angord,scatrecoil,allhad,
compress
5559 DOUBLE PRECISION z,
q
5572 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
5575 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
5576 LOGICAL angord,scatrecoil,allhad,
compress
5578 DOUBLE PRECISION z,
q
5591 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
5594 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
5595 LOGICAL angord,scatrecoil,allhad,
compress
5597 DOUBLE PRECISION z,
q,
ei
5600 & - 2.*lps**4*
ei(2.*(
log(
q**2/lps**2)+
log(
z)))/
q**4
5601 & + 2.*lps**6*
ei(3.*(
log(
q**2/lps**2)+
log(
z)))/
q**6)/
5613 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
5616 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
5617 LOGICAL angord,scatrecoil,allhad,
compress
5619 DOUBLE PRECISION z,
q,
ei
5622 & - 2.*lps**4*
ei(2.*(
log(
q**2/lps**2)+
log(1.-
z)))/
q**4
5623 & + 2.*lps**6*
ei(3.*(
log(
q**2/lps**2)+
log(1.-
z)))/
q**6)/
5632 DOUBLE PRECISION FUNCTION gett(MINT,MAXT,MDEB)
5635 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
5638 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
5639 LOGICAL angord,scatrecoil,allhad,
compress
5641 DOUBLE PRECISION tmin,tmax,maxi,
pyr,
r1,
r2,
alphas,
pi,
y,maxt,
5643 DATA pi/3.141592653589793d0/
5647 IF(tmin.GT.tmax)
THEN
5652 t=tmax*tmin/(tmax+
r1*(tmin-tmax))
5667 DOUBLE PRECISION FUNCTION ei(X)
5670 common/hepmcid/hpmcfid,logfid
5671 integer hpmcfid,logfid
5673 common/expint/eixs(3,1000),valmax,nval
5675 DOUBLE PRECISION eixs,valmax
5678 DOUBLE PRECISION x,
r,ga,xa(2),ya(2),
y,
dy,a,b
5679 DOUBLE PRECISION ystart,epsi,hfirst
5682 IF(dabs(
x).GT.valmax)
5683 &
write(logfid,*)
'warning: value out of array in Ei(x)',
x,valmax
5686 lmax=int(
x*nval/valmax)
5689 IF((
line.GT.999).OR.(
line.LT.1))
THEN
5690 write(logfid,*)
'ERROR in EI: line number out of bound',
line
5696 a=(ya(2)-ya(1))/(xa(2)-xa(1))
5700 lmax=int(-
x*nval/valmax)
5703 IF((
line.GT.999).OR.(
line.LT.1))
THEN
5704 write(logfid,*)
'ERROR in EI: line number out of bound',
line
5710 a=(ya(2)-ya(1))/(xa(2)-xa(1))
5722 DOUBLE PRECISION FUNCTION pqq(Z)
5725 pqq=4.*(1.+
z**2)/(3.*(1.-
z))
5733 DOUBLE PRECISION FUNCTION pgq(Z)
5736 pgq=4.*(1.+(1.-
z)**2)/(3.*
z)
5744 DOUBLE PRECISION FUNCTION pgg(Z)
5755 DOUBLE PRECISION FUNCTION pqg(Z)
5758 pqg=0.5*(
z**2 + (1.-
z)**2)
5769 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
5772 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
5773 LOGICAL angord,scatrecoil,allhad,
compress
5775 DOUBLE PRECISION t,l0,
pi,lambda
5776 DATA pi/3.141592653589793d0/
5789 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
5792 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
5793 LOGICAL angord,scatrecoil,allhad,
compress
5795 common/splitint/splitiggv(1000,1000),splitiqqv(1000,1000),
5796 &splitiqgv(1000,1000),qval(1000),zmval(1000),qmax,zmmin,npoint
5798 DOUBLE PRECISION splitiggv,splitiqqv,splitiqgv,
5799 &qval,zmval,qmax,zmmin
5801 common/intsplitf/qquad,fm
5802 DOUBLE PRECISION qquad,fm
5808 DOUBLE PRECISION emax,zmmax,epsi,hfirst,ystart,lnzmmin,
5809 &lnzmmax,zm,zm2,
q,
getmsmax,avmom(5),shat,pcms2
5814 call
maxscatcen(avmom(1),avmom(2),avmom(3),avmom(4),avmom(5))
5815 shat = avmom(5)**2 +
5816 & 2.*emax*(avmom(4)+sqrt(avmom(1)**2+avmom(2)**2+avmom(3)**2))
5817 pcms2 = (shat-avmom(5)**2)**2/(4.*shat)
5818 qmax = sqrt(scalefacm*4.*pcms2)
5828 q=(
i-1)*(qmax-q0)/nstep+q0
5832 zm=exp((
j-1)*(lnzmmax-lnzmmin)/nstep+lnzmmin)
5834 IF(
q**2.LT.q0**2)
THEN
5837 zm2=0.5-0.5*sqrt(1.-q0**2/
q**2)
5848 CALL
odeint(ystart,zm,1.-zm,epsi,hfirst,0d0,2)
5849 splitiqqv(
i,
j)=ystart
5853 CALL
odeint(ystart,zm,1.-zm,epsi,hfirst,0d0,3)
5854 splitiggv(
i,
j)=ystart
5858 CALL
odeint(ystart,zm,1.-zm,epsi,hfirst,0d0,4)
5859 splitiqgv(
i,
j)=ystart
5874 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
5877 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
5878 LOGICAL angord,scatrecoil,allhad,
compress
5880 common/pdfs/qinqx(2,1000),ginqx(2,1000),qingx(2,1000),
5882 DOUBLE PRECISION qinqx,ginqx,qingx,gingx
5884 common/pdfintv/
xmax,
z
5885 DOUBLE PRECISION xmax,
z
5892 &q2max,deltaq2,avmom(5),shat,pcms2
5895 call
maxscatcen(avmom(1),avmom(2),avmom(3),avmom(4),avmom(5))
5896 shat = avmom(5)**2 +
5897 & 2.*emax*(avmom(4)+sqrt(avmom(1)**2+avmom(2)**2+avmom(3)**2))
5898 pcms2 = (shat-avmom(5)**2)**2/(4.*shat)
5899 q2max = scalefacm*4.*pcms2
5901 deltaq2=
log(q2max)-
log(q0**2)
5911 q2 = exp((
j-1)*deltaq2/999.d0 +
log(q0**2))
5931 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
5934 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
5935 LOGICAL angord,scatrecoil,allhad,
compress
5937 common/xsecs/intq1(1001,101),intq2(1001,101),
5938 &intg1(1001,101),intg2(1001,101)
5939 DOUBLE PRECISION intq1,intq2,intg1,intg2
5941 common/xsecv/qlow,mdx
5942 DOUBLE PRECISION qlow,mdx
5948 DOUBLE PRECISION emax,tmax,tmaxmax,deltatmax,ystart,hfirst,epsi,
5952 call
maxscatcen(avmom(1),avmom(2),avmom(3),avmom(4),avmom(5))
5953 shat = avmom(5)**2 +
5954 & 2.*emax*(avmom(4)+sqrt(avmom(1)**2+avmom(2)**2+avmom(3)**2))
5955 pcms2 = (shat-avmom(5)**2)**2/(4.*shat)
5956 tmaxmax = scalefacm*4.*pcms2
5957 deltatmax=(
log(tmaxmax)-
5958 &
log(q0**2*(1.d0+1.d-6)/scalefacm**2))/999.d0
5961 deltamd=(mdmax-mdmin)/99.d0
5964 tmax = exp((
j-1)*deltatmax
5965 & +
log(q0**2*(1.d0+1.d-6)/scalefacm**2))
5971 mdx=mdmin+(
k-1)*deltamd
5976 IF(tmax.LT.q0**2/scalefacm**2)
THEN
5984 hfirst=0.01*(tmax-q0**2/scalefacm**2)
5986 CALL
odeint(ystart,q0**2/scalefacm**2,tmax,epsi,hfirst
5991 hfirst=0.01*(tmax-q0**2/scalefacm**2)
5993 CALL
odeint(ystart,q0**2/scalefacm**2,tmax,epsi,hfirst
5999 CALL
odeint(ystart,q0**2/scalefacm**2,tmax,epsi,hfirst
6005 CALL
odeint(ystart,q0**2/scalefacm**2,tmax,epsi,hfirst
6021 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
6024 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
6025 LOGICAL angord,scatrecoil,allhad,
compress
6027 common/insuda/sudaqq(1000,2),sudaqg(1000,2),sudagg(1000,2),
6029 DOUBLE PRECISION sudaqq,sudaqg,sudagg,sudagc
6038 call
maxscatcen(avmom(1),avmom(2),avmom(3),avmom(4),avmom(5))
6039 shat = avmom(5)**2 +
6040 & 2.*emax*(avmom(4)+sqrt(avmom(1)**2+avmom(2)**2+avmom(3)**2))
6041 pcms2 = (shat-avmom(5)**2)**2/(4.*shat)
6042 qmax = sqrt(scalefacm*4.*pcms2)
6043 delta=(
log(3.*qmax)-
log(q0**2*(1.d0+1.d-6)))/999.d0
6045 q = exp((
i-1)*delta +
log(q0**2*(1.d0+1.d-6)))
6065 common/expint/eixs(3,1000),valmax,nval
6067 DOUBLE PRECISION eixs,valmax
6070 DOUBLE PRECISION x,epsi,hfirst,ystart,
ei,ga,
r
6077 x=
i*valmax/(nval*1.d0)
6082 CALL
odeint(ystart,
x,1000.d0,epsi,hfirst,0.d0,5)
6087 ELSE IF (
x.LE.40.0)
THEN
6091 r=
r*
k*
x/(
k+1.0d0)**2
6093 IF (dabs(
r/
ei).LE.1.0d-15)
go to 20
6095 20 ga=0.5772156649015328d0
6117 common/hepmcid/hpmcfid,logfid
6118 integer hpmcfid,logfid
6120 integer nmax,nstep,w1
6121 double precision ystart,a,b,
eps,
h1,hmin,
x,
h,
y,dydx,
6122 &
deriv,yscale,hdid,hnew
6130 yscale = abs(
y) + abs(
h*dydx) + 1.
e-25
6131 if (((
x +
h - b)*
h).gt.0.)
h = b-
x
6133 if ((
x - b)*
h.ge.0)
then
6138 if (abs(
h).lt.abs(hmin))
then
6139 write(logfid,*)
'Error in odeint: stepsize too small',w1
6144 write(logfid,*)
'Error in odeint: too many steps',w1
6153 subroutine rkstepper(x,y,dydx,htest,hdid,hnew,yscale,eps,w1)
6156 common/hepmcid/hpmcfid,logfid
6157 integer hpmcfid,logfid
6160 double precision x,
y,dydx,htest,hdid,hnew,yscale,
eps,
6161 &yhalf,y1,y2,
rk4step,dydxhalf,xnew,delta,err,
h,safety, powerdown,
6162 &powerup,maxup,maxdown,
deriv,fac
6164 data powerdown/0.25/
6174 write(logfid,*)
'Error in rkstepper: step size not significant'
6182 err = abs(delta)/(yscale*
eps)
6185 fac =
max(1./maxdown,safety/err**powerdown)
6192 fac = min(maxup,safety/err**powerup)
6225 common/hepmcid/hpmcfid,logfid
6226 integer hpmcfid,logfid
6230 DOUBLE PRECISION p,
v
6232 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
6235 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
6236 LOGICAL angord,scatrecoil,allhad,
compress
6238 common/time/mv(23000,5)
6247 INTEGER line,
i,nnull
6248 DOUBLE PRECISION dtmax,sigmamax,neffmax,linvmax,
pyr,
6249 &
r,toff,xs,ys,zs,
ts,
getsscat,
getmsmax,
getmdmin,msmax,mdmin,
6252 &sigmamin,neffmin,tstart,dtmax1,deltat
6260 IF(
k(
line,2).EQ.21)
THEN
6270 IF((tstart+dtmax).GE.ltime)dtmax=ltime-tstart
6271 IF(dtmax.LT.0.d0)
RETURN
6274 toff=tstart-mv(
line,4)
6282 &
p(
line,5),0.d0,ptype,
'C',xstart,ystart,zstart,tstart,1)
6285 &
p(
line,5),0.d0,ptype,
'C',xstart,ystart,zstart,tstart,2)
6288 linvmax=5.d0*
max(neffmin*sigmamax,neffmax*sigmamin)
6289 if(linvmax.eq.0.d0)
return
6292 deltat=deltat-
log(
pyr(0))/linvmax
6300 tau=sqrt(
ts**2-zs**2)
6303 IF((
tau.GT.1.d0).AND.(neff.EQ.0.d0))
THEN
6312 IF((deltat.GT.dtmax).OR.stopnow)
THEN
6316 IF(neff.GT.0.d0)
THEN
6318 &
p(
line,5),0.d0,ptype,
'C',xs,ys,zs,
ts,0)
6323 IF(
weight.GT.1.d0+1d-6)
then
6324 if (
line.ne.errl)
then
6325 write(logfid,*)
'error in GETDELTAT: weight > 1',
weight,
6326 & neff*
sigma/(neffmax*sigmamin),neff*
sigma/(neffmin*sigmamax),
6343 double precision lambda,disc,
p,
pyr,u,
v,
pi
6344 data pi/3.141592653589793d0/
6346 if (lambda.gt.745.d0)
then
6350 & int(sqrt(lambda)*sqrt(-2.*
log(u))*cos(2.*
pi*
v)+lambda)
6372 IF(abs(
id).LT.100)
THEN
6375 IF(mod(int(abs(
id)/10.),10).EQ.0)
THEN
6392 IF(abs(
id).LT.1000)
THEN
6395 IF(mod(int(
id/10),10).EQ.0)
THEN
6410 IF((abs(
id).EQ.11).OR.(abs(
id).EQ.13).OR.(abs(
id).EQ.15))
THEN
6441 DOUBLE PRECISION p,
v
6445 if ((
k(l,2).ne.91).and.(
k(l,2).ne.92))
then
6449 if ((
k(
k(l,3),3).eq.0).or.(
isparton(
k(
k(
k(l,3),3),2))))
then
6465 DOUBLE PRECISION p,
v
6469 if ((
k(l,2).ne.91).and.(
k(l,2).ne.92))
then
6493 DOUBLE PRECISION p,
v
6497 if (((
k(
k(l,3),2).EQ.91).OR.(
k(
k(l,3),2).EQ.92))
6515 common/hepmcid/hpmcfid,logfid
6516 integer hpmcfid,logfid
6519 DOUBLE PRECISION p,
v
6521 common/angor/za(23000),zd(23000),thetaa(23000),qqbard(23000)
6522 DOUBLE PRECISION za,zd,thetaa
6525 common/time/mv(23000,5)
6528 common/colour/trip(23000),anti(23000),colmax
6529 INTEGER trip,anti,colmax
6531 integer l1,
i,
j,nold,nnew,nstart
6536 if (((
k(
i,1).eq.11).or.(
k(
i,1).eq.12).or.(
k(
i,1).eq.13)).and.
6546 do 779
i=nstart,nold
6547 if (((
k(
i,1).ne.11).and.(
k(
i,1).ne.12).and.(
k(
i,1).ne.13)).or.
6558 thetaa(nnew)=thetaa(
i)
6559 qqbard(nnew)=qqbard(
i)
6565 if (l1.eq.
i) l1=nnew
6570 if ((nold-
n).le.10)
then
6589 if (
n.gt.23000)
write(logfid,*)
'Error in compressevent: n = ',
n
6590 if (l1.gt.
n)
write(logfid,*)
'Error in compressevent: l1 = ',l1
6603 common/hepmcid/hpmcfid,logfid
6604 integer hpmcfid,logfid
6607 DOUBLE PRECISION p,
v
6609 common/angor/za(23000),zd(23000),thetaa(23000),qqbard(23000)
6610 DOUBLE PRECISION za,zd,thetaa
6613 common/time/mv(23000,5)
6616 common/colour/trip(23000),anti(23000),colmax
6617 INTEGER trip,anti,colmax
6627 IF(
col)
write(logfid,*)
i,
' (',trip(
i),
',',anti(
i),
') [',
6628 &
k(
i,3),
k(
i,4),
k(
i,5),
' ] {',
k(
i,2),
k(
i,1),
' } ',
6644 DOUBLE PRECISION p,
v
6649 common/
param/q0,lps,lqcd,ltime,scalefacm,angord,scatrecoil,
6652 DOUBLE PRECISION q0,lqcd,ltime,lps,scalefacm
6653 LOGICAL angord,scatrecoil,allhad,
compress
6655 common/evrecord/nsim,
npart,
offset,hadrotype,sqrts,collider,hadro,
6656 &shorthepmc,
channel,isochannel
6658 double precision sqrts
6660 character*2 isochannel
6661 logical hadro,shorthepmc
6663 common/storescatcen/nscatcen,maxnscatcen,scatflav(10000),
6664 &scatcen(10000,5),writescatcen,writedummies
6665 integer nscatcen,maxnscatcen,scatflav
6666 double precision scatcen
6667 logical writescatcen,writedummies
6669 INTEGER evnum,pbarcode,vbarcode,codelist(25000),
i,pid,nstart,
6670 &nfirst,nvertex,ntot,
j,codefirst
6671 DOUBLE PRECISION mproton,mneutron,pdummy,pscatcen
6674 character*2 beam1,beam2
6675 data mproton/0.9383/
6676 data mneutron/0.9396/
6679 5000
FORMAT(a2,i10,
i3,3e14.6,2
i2,i6,4
i2,e14.6)
6680 5100
FORMAT(a2,2e14.6)
6681 5200
FORMAT(a2,6i7,2
i2,1i7,4e14.6)
6682 5300
FORMAT(a2,2
i2,5e14.6,2
i2)
6683 5400
FORMAT(a2,i6,6
i2,i6,
i2)
6684 5500
FORMAT(a2,i6,i6,5e14.6,3
i2,i6,
i2)
6689 if (shorthepmc)
then
6691 IF(collider.EQ.
'EEJJ')
THEN
6700 if (((
k(
i,1).lt.6).or.(
k(
i,1).eq.17)))
6703 if(writescatcen) nfirst=nfirst+nscatcen
6704 if(writedummies) nfirst=nfirst+nscatcen
6706 WRITE(
j,5000)
'E ',evnum,-1,0.d0,0.d0,0.d0,0,0,nvertex,1,2,0,1,
6708 WRITE(
j,
'(A2,I2,A5)')
'N ',1,
'"0"'
6709 WRITE(
j,
'(A)')
'U GEV MM'
6710 WRITE(
j,5100)
'C ',
pari(1)*1.d9,0.d0
6711 WRITE(
j,5200)
'H ',0,0,0,0,0,0,0,0,0,0.d0,0.d0,0.d0,0.d0
6712 WRITE(
j,5300)
'F ',0,0,-1.d0,-1.d0,-1.d0,-1.d0,-1.d0,0,0
6714 IF(collider.EQ.
'EEJJ')
THEN
6715 WRITE(
j,5400)
'V ',-1,0,0,0,0,0,2,1,0
6716 WRITE(
j,5500)
'P ',1,-11,0.d0,0.d0,sqrts/2.,sqrts/2.,
6717 & 0.00051,2,0,0,-1,0
6718 WRITE(
j,5500)
'P ',2,11,0.d0,0.d0,-sqrts/2.,sqrts/2.,
6719 & 0.00051,2,0,0,-1,0
6720 WRITE(
j,5500)
'P ',3,23,0.d0,0.d0,0.d0,sqrts,
6722 WRITE(
j,5400)
'V ',-2,0,0,0,0,0,0,2,0
6723 WRITE(
j,5500)
'P ',4,pid,sqrts/2.,0.d0,0.d0,sqrts/2.,
6725 WRITE(
j,5500)
'P ',5,-pid,-sqrts/2.,0.d0,0.d0,sqrts/2.,
6727 WRITE(
j,5400)
'V ',-3,0,0,0,0,0,0,nfirst,0
6729 WRITE(
j,5400)
'V ',-1,0,0,0,0,0,2,nfirst,0
6730 if (beam1.eq.
'p+')
then
6731 WRITE(
j,5500)
'P ',1,2212,0.d0,0.d0,
6732 & sqrt(sqrts**2/4.-mproton**2),sqrts/2.,mproton,2,0,0,-1,0
6734 WRITE(
j,5500)
'P ',1,2112,0.d0,0.d0,
6735 & sqrt(sqrts**2/4.-mneutron**2),sqrts/2.,mneutron,2,0,0,-1,0
6737 if (beam2.eq.
'p+')
then
6738 WRITE(
j,5500)
'P ',2,2212,0.d0,0.d0,
6739 & -sqrt(sqrts**2/4.-mproton**2),sqrts/2.,mproton,2,0,0,-1,0
6741 WRITE(
j,5500)
'P ',2,2112,0.d0,0.d0,
6742 & -sqrt(sqrts**2/4.-mneutron**2),sqrts/2.,mneutron,2,0,0,-1,0
6746 if(writescatcen)
then
6749 WRITE(
j,5500)
'P ',pbarcode,scatflav(
i),scatcen(
i,1),
6750 & scatcen(
i,2),scatcen(
i,3),scatcen(
i,4),scatcen(
i,5),
6755 if(writedummies)
then
6758 pscatcen=sqrt(scatcen(
i,1)**2+scatcen(
i,2)**2+
6760 WRITE(
j,5500)
'P ',pbarcode,111,pdummy*scatcen(
i,1)/pscatcen,
6761 & pdummy*scatcen(
i,2)/pscatcen,pdummy*scatcen(
i,3)/pscatcen,
6762 & pdummy,0.d0,1,0,0,0,0
6767 if(((
k(
i,1).lt.6).or.(
k(
i,1).eq.17)))
then
6769 if((
k(
i,1).eq.3).or.(
k(
i,1).eq.5))
then
6770 WRITE(
j,5500)
'P ',pbarcode,
k(
i,2),
p(
i,1),
p(
i,2),
p(
i,3),
6771 &
p(
i,4),
p(
i,5),4,0,0,0,0
6773 WRITE(
j,5500)
'P ',pbarcode,
k(
i,2),
p(
i,1),
p(
i,2),
p(
i,3),
6774 &
p(
i,4),
p(
i,5),1,0,0,0,0
6784 IF(collider.EQ.
'EEJJ')
THEN
6801 & .AND.(
k(
i,4).NE.0)) nvertex=nvertex+1
6805 if(writescatcen) nfirst=nfirst+nscatcen
6806 if(writedummies) nfirst=nfirst+nscatcen
6808 WRITE(
j,5000)
'E ',evnum,-1,0.d0,0.d0,0.d0,0,0,nvertex,
6810 WRITE(
j,
'(A2,I2,A5)')
'N ',1,
'"0"'
6811 WRITE(
j,
'(A)')
'U GEV MM'
6812 WRITE(
j,5100)
'C ',
pari(1)*1.d9,0.d0
6813 WRITE(
j,5200)
'H ',0,0,0,0,0,0,0,0,0,0.d0,0.d0,0.d0,0.d0
6814 WRITE(
j,5300)
'F ',0,0,-1.d0,-1.d0,-1.d0,-1.d0,-1.d0,0,0
6817 IF(collider.EQ.
'EEJJ')
THEN
6824 IF(collider.EQ.
'EEJJ')
THEN
6825 WRITE(
j,5400)
'V ',-1,0,0,0,0,0,2,1,0
6826 WRITE(
j,5500)
'P ',1,-11,0.d0,0.d0,sqrts/2.,sqrts/2.,
6827 & 0.00051,2,0,0,-1,0
6828 WRITE(
j,5500)
'P ',2,11,0.d0,0.d0,-sqrts/2.,sqrts/2.,
6829 & 0.00051,2,0,0,-1,0
6830 WRITE(
j,5500)
'P ',3,23,0.d0,0.d0,0.d0,sqrts,
6832 WRITE(
j,5400)
'V ',-2,0,0,0,0,0,0,2,0
6833 WRITE(
j,5500)
'P ',4,pid,sqrts/2.,0.d0,0.d0,sqrts/2.,
6835 WRITE(
j,5500)
'P ',5,-pid,-sqrts/2.,0.d0,0.d0,sqrts/2.,
6837 WRITE(
j,5400)
'V ',vbarcode,0,0,0,0,0,0,nfirst,0
6839 WRITE(
j,5400)
'V ',-1,0,0,0,0,0,2,nfirst,0
6840 if (beam1.eq.
'p+')
then
6841 WRITE(
j,5500)
'P ',1,2212,0.d0,0.d0,
6842 & sqrt(sqrts**2/4.-mproton**2),sqrts/2.,mproton,2,0,0,-1,0
6844 WRITE(
j,5500)
'P ',1,2112,0.d0,0.d0,
6845 & sqrt(sqrts**2/4.-mneutron**2),sqrts/2.,mneutron,2,0,0,-1,0
6847 if (beam2.eq.
'p+')
then
6848 WRITE(
j,5500)
'P ',2,2212,0.d0,0.d0,
6849 & -sqrt(sqrts**2/4.-mproton**2),sqrts/2.,mproton,2,0,0,-1,0
6851 WRITE(
j,5500)
'P ',2,2112,0.d0,0.d0,
6852 & -sqrt(sqrts**2/4.-mneutron**2),sqrts/2.,mneutron,2,0,0,-1,0
6856 codefirst=nfirst+pbarcode
6859 if(writescatcen)
then
6862 WRITE(
j,5500)
'P ',pbarcode,scatflav(
i),scatcen(
i,1),
6863 & scatcen(
i,2),scatcen(
i,3),scatcen(
i,4),scatcen(
i,5),
6868 if(writedummies)
then
6871 pscatcen=sqrt(scatcen(
i,1)**2+scatcen(
i,2)**2+
6873 WRITE(
j,5500)
'P ',pbarcode,111,pdummy*scatcen(
i,1)/pscatcen,
6874 & pdummy*scatcen(
i,2)/pscatcen,pdummy*scatcen(
i,3)/pscatcen,
6875 & pdummy,0.d0,1,0,0,0,0
6884 IF (pbarcode.EQ.codefirst) goto 130
6889 codelist(
i)=vbarcode
6890 WRITE(
j,5500)
'P ',pbarcode,
k(
i,2),
p(
i,1),
p(
i,2),
p(
i,3),
6891 &
p(
i,4),
p(
i,5),2,0,0,vbarcode,0
6893 WRITE(
j,5500)
'P ',pbarcode,
k(
i,2),
p(
i,1),
p(
i,2),
p(
i,3),
6894 &
p(
i,4),
p(
i,5),1,0,0,0,0
6907 codelist(
i)=codelist(
k(
i,3))
6911 IF((
k(
i,3).NE.
k(
i-1,3)))
THEN
6913 WRITE(
j,5400)
'V ',codelist(
k(
i,3)),0,0,0,0,0,0,
6914 &
k(
k(
i,3),5)-
k(
k(
i,3),4)+1,0
6919 codelist(
i)=vbarcode
6920 WRITE(
j,5500)
'P ',pbarcode,
k(
i,2),
p(
i,1),
p(
i,2),
p(
i,3),
6921 &
p(
i,4),
p(
i,5),2,0,0,vbarcode,0
6923 WRITE(
j,5500)
'P ',pbarcode,
k(
i,2),
p(
i,1),
p(
i,2),
p(
i,3),
6924 &
p(
i,4),
p(
i,5),1,0,0,0,0
6945 write(
fid,*)
' _______________'//
6946 &
'__________________________ '
6949 write(
fid,*)
' | JJJJJ EEEEE '//
6951 write(
fid,*)
' | J E '//
6953 write(
fid,*)
' _________________| J EEE '//
6954 &
' W W W EEE L |_________________ '
6955 write(
fid,*)
'| | J J E '//
6957 write(
fid,*)
'| | JJJ EEEEE '//
6958 &
' W W EEEEE LLLLL | |'
6959 write(
fid,*)
'| |_______________'//
6960 &
'__________________________| |'
6964 &
'this is JEWEL 2.1.0 |'
6967 write(
fid,*)
'| Copyright Korinna C. Zapp (2016)'//
6968 &
' [Korinna.Zapp@cern.ch] |'
6971 write(
fid,*)
'| The JEWEL homepage is jewel.hepforge.org '//
6975 write(
fid,*)
'| The medium model was partly '//
6976 &
'implemented by Jochen Klein |'
6977 write(
fid,*)
'| [Jochen.Klein@cern.ch]. Raghav '//
6978 &
'Kunnawalkam Elayavalli helped with the |'
6979 write(
fid,*)
'| implementation of the V+jet processes '//
6980 &
'[raghav.k.e@cern.ch]. |'
6983 write(
fid,*)
'| Please cite JHEP 1303 (2013) '//
6984 &
'080 [arXiv:1212.1599] and optionally |'
6985 write(
fid,*)
'| EPJC C60 (2009) 617 [arXiv:0804.3568] '//
6986 &
'for the physics and arXiv:1311.0048 |'
6987 write(
fid,*)
'| for the code. The reference for '//
6988 &
'V+jet processes is EPJC 76 (2016) no.12 695 |'
6989 write(
fid,*)
'| [arXiv:1608.03099] and for recoil effects'//
6990 &
' it is arXiv:1707.01539. |'
6993 write(
fid,*)
'| JEWEL contains code provided by '//
6994 &
'S. Zhang and J. M. Jing |'
6995 write(
fid,*)
'| (Computation of Special Functions, '//
6996 &
'John Wiley & Sons, New York, 1996 and |'
6997 write(
fid,*)
'| http://jin.ece.illinois.edu) for '//
6998 &
'computing the exponential integral Ei(x). |'
7001 write(
fid,*)
'| JEWEL relies heavily on PYTHIA 6'//
7002 &
' for the event generation. The modified |'
7003 write(
fid,*)
'| version of PYTHIA 6.4.25 that is'//
7004 &
' shipped with JEWEL is, however, not an |'
7005 write(
fid,*)
'| official PYTHIA release and must'//
7006 &
' not be used for anything else. Please |'
7007 write(
fid,*)
'| refer to results as "JEWEL+PYTHIA".'//
7011 write(
fid,*)
'|_________________________________'//
7012 &
'____________________________________________|'
7024 common/hepmcid/hpmcfid,logfid
7025 integer hpmcfid,logfid
7027 integer*4 date(3),time(3)
7029 1000
format (
i2.2,
'.',
i2.2,
'.',
i4.4,
', ',
7030 &
i2.2,
':',
i2.2,
':',
i2.2 )
7033 write(logfid,1000)date,time