Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F7877351
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
198 KB
Subscribers
None
View Options
Index: branches/fragphotons/jewel-2.2.0.f
===================================================================
--- branches/fragphotons/jewel-2.2.0.f (revision 390)
+++ branches/fragphotons/jewel-2.2.0.f (revision 391)
@@ -1,6942 +1,6964 @@
PROGRAM JEWEL
IMPLICIT NONE
C--Common block of Pythia
COMMON/PYJETS/N,NPAD,K(23000,5),P(23000,5),V(23000,5)
INTEGER N,NPAD,K
DOUBLE PRECISION P,V
COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
INTEGER MSTU,MSTJ
DOUBLE PRECISION PARU,PARJ
COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
INTEGER MDCY,MDME,KFDP
DOUBLE PRECISION BRAT
COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
INTEGER MSEL,MSELPD,MSUB,KFIN
DOUBLE PRECISION CKIN
COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
INTEGER MSTP,MSTI
DOUBLE PRECISION PARP,PARI
COMMON/PYDATR/MRPY(6),RRPY(100)
INTEGER MRPY
DOUBLE PRECISION RRPY
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
C--use nuclear pdf?
COMMON/NPDF/MASS,NSET,EPS09,INITSTR
INTEGER NSET
DOUBLE PRECISION MASS
LOGICAL EPS09
CHARACTER*10 INITSTR
C--number of protons
common/np/nproton
integer nproton
C--organisation of event record
common/evrecord/nsim,npart,offset,hadrotype,sqrts,collider,hadro,
&shorthepmc,channel,isochannel
integer nsim,npart,offset,hadrotype
double precision sqrts
character*4 collider,channel
character*2 isochannel
logical hadro,shorthepmc
C--discard event flag
COMMON/DISC/NDISC,NSTRANGE,NGOOD,errcount,wdisc,DISCARD
LOGICAL DISCARD
INTEGER NDISC,NSTRANGE,NGOOD,errcount
double precision wdisc
C--event weight
COMMON/WEIGHT/EVWEIGHT,sumofweights
double precision EVWEIGHT,sumofweights
C--number of scattering events
COMMON/CHECK/NSCAT,NSCATEFF,NSPLIT,nspliti,nsplitf,nistry,
&nisfail,nfsfail,nfstry,nttot,ntrej
DOUBLE PRECISION NSCAT,NSCATEFF,NSPLIT,nspliti,nsplitf,nistry,
&nisfail,nfsfail,nfstry,nttot,ntrej
C--number of extrapolations in tables
common/extrapolations/ntotspliti,noverspliti,ntotpdf,noverpdf,
&ntotxsec,noverxsec,ntotsuda,noversuda
integer ntotspliti,noverspliti,ntotpdf,noverpdf,
&ntotxsec,noverxsec,ntotsuda,noversuda
C--local variables
integer j,i,kk,poissonian
integer nsimpp,nsimpn,nsimnp,nsimnn,nsimsum,nsimchn
double precision sumofweightstot,wdisctot,scalefac
double precision gettemp,r,tau
character*2 b1,b2
call init()
SUMOFWEIGHTSTOT=0.d0
WDISCTOT=0.d0
C--e+ + e- event generation
if (collider.eq.'EEJJ') then
b1 = 'e+'
b2 = 'e-'
write(logfid,*)
write(logfid,*)
&'####################################################'
write(logfid,*)
write(logfid,*)'generating ',nsim,' events in ',b1,' + ',b2,
&' channel'
write(logfid,*)
write(logfid,*)
&'####################################################'
write(logfid,*)
SUMOFWEIGHTS=0.d0
WDISC=0.d0
call initpythia(b1,b2)
write(logfid,*)
C--e+ + e- event loop
DO 100 J=1,NSIM
call genevent(j,b1,b2)
100 CONTINUE
sumofweightstot = sumofweightstot+sumofweights
wdisctot = wdisctot + wdisc
write(logfid,*)
write(logfid,*)'cross section in e+ + e- channel:',PARI(1),'mb'
write(logfid,*)'sum of event weights in e+ + e- channel:',
& sumofweights-wdisc
write(logfid,*)
else
C--hadronic event generation
if (isochannel.eq.'PP') then
nsimpp = nsim
nsimpn = 0
nsimnp = 0
nsimnn = 0
elseif (isochannel.eq.'PN') then
nsimpp = 0
nsimpn = nsim
nsimnp = 0
nsimnn = 0
elseif (isochannel.eq.'NP') then
nsimpp = 0
nsimpn = 0
nsimnp = nsim
nsimnn = 0
elseif (isochannel.eq.'NN') then
nsimpp = 0
nsimpn = 0
nsimnp = 0
nsimnn = nsim
else
nsimpp = poissonian(nsim*nproton**2/mass**2)
nsimpn = poissonian(nsim*nproton*(mass-nproton*1.d0)/mass**2)
nsimnp = poissonian(nsim*nproton*(mass-nproton*1.d0)/mass**2)
nsimnn = poissonian(nsim*(mass-nproton*1.d0)**2/mass**2)
nsimsum = nsimpp + nsimpn + nsimnp + nsimnn
scalefac = nsim*1.d0/(nsimsum*1.d0)
nsimpp = int(nsimpp*scalefac)
nsimpn = int(nsimpn*scalefac)
nsimnp = int(nsimnp*scalefac)
nsimnn = int(nsimnn*scalefac)
nsimsum = nsimpp + nsimpn + nsimnp + nsimnn
endif
C--loop over channels
do 101 kk=1,4
if (kk.eq.1) then
b1 = 'p+'
b2 = 'p+'
nsimchn = nsimpp
elseif (kk.eq.2) then
b1 = 'p+'
b2 = 'n0'
nsimchn = nsimpn
elseif (kk.eq.3) then
b1 = 'n0'
b2 = 'p+'
nsimchn = nsimnp
else
b1 = 'n0'
b2 = 'n0'
nsimchn = nsimnn
endif
write(logfid,*)
write(logfid,*)
&'####################################################'
write(logfid,*)
write(logfid,*)'generating ',nsimchn,' events in ',
&b1,' + ',b2,' channel'
write(logfid,*)
write(logfid,*)
&'####################################################'
write(logfid,*)
SUMOFWEIGHTS=0.d0
WDISC=0.d0
call initpythia(b1,b2)
write(logfid,*)
C--event loop
DO 102 J=1,nsimchn
call genevent(j,b1,b2)
102 CONTINUE
sumofweightstot = sumofweightstot+sumofweights
wdisctot = wdisctot + wdisc
write(logfid,*)
write(logfid,*)'cross section in ',b1,' + ',b2,' channel:',
& PARI(1),'mb'
write(logfid,*)'sum of event weights in ',b1,' + ',b2,
& ' channel:',sumofweights-wdisc
write(logfid,*)
101 continue
endif
C--finish
WRITE(HPMCFID,'(A)')'HepMC::IO_GenEvent-END_EVENT_LISTING'
WRITE(HPMCFID,*)
CLOSE(HPMCFID,status='keep')
write(logfid,*)
write(logfid,*)'mean number of scatterings:',
& NSCAT/(SUMOFWEIGHTSTOT-WDISCTOT)
write(logfid,*)'mean number of effective scatterings:',
& NSCATEFF/(SUMOFWEIGHTSTOT-WDISCTOT)
write(logfid,*)'mean number of splittings:',
& NSPLIT/(SUMOFWEIGHTSTOT-WDISCTOT)
write(logfid,*)'mean number of splittings from IS shower:',
& nspliti/(SUMOFWEIGHTSTOT-WDISCTOT)
write(logfid,*)'mean number of splittings from FS shower:',
& nsplitf/(SUMOFWEIGHTSTOT-WDISCTOT)
write(logfid,*)'fraction of rejected IS splittings:',
& nisfail/nistry
write(logfid,*)'fraction of rejected FS splittings:',
& nfsfail/nfstry
write(logfid,*)'fraction of rejected momentum transfers:',
& ntrej/nttot
write(logfid,*)
write(logfid,*)'number of extrapolations in splitting integral: ',
& noverspliti,' (',(noverspliti*1.d0)/(ntotspliti*1.d0),'%)'
write(logfid,*)
& 'number of extrapolations in splitting partonic PDFs: ',
& noverpdf,' (',(noverpdf*1.d0)/(ntotpdf*1.d0),'%)'
write(logfid,*)
& 'number of extrapolations in splitting cross sections: ',
& noverxsec,' (',(noverxsec*1.d0)/(ntotxsec*1.d0),'%)'
write(logfid,*)
& 'number of extrapolations in Sudakov form factor: ',
& noversuda,' (',(noversuda*1.d0)/(ntotsuda*1.d0),'%)'
write(logfid,*)
write(logfid,*)'number of good events: ',ngood
write(logfid,*)'total number of discarded events: ',NDISC
write(logfid,*)'number of events for which conversion '//
&'to hepmc failed: ',NSTRANGE
call printtime
close(logfid,status='keep')
END
***********************************************************************
***********************************************************************
*** END OF MAIN PROGRAM - NOW COME THE SUBROUTINES ****************
***********************************************************************
***********************************************************************
***********************************************************************
*** subroutine init
***********************************************************************
subroutine init()
implicit none
INTEGER PYCOMP
INTEGER NMXHEP
C--Common block of Pythia
COMMON/PYJETS/N,NPAD,K(23000,5),P(23000,5),V(23000,5)
INTEGER N,NPAD,K
DOUBLE PRECISION P,V
COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
INTEGER MSTU,MSTJ
DOUBLE PRECISION PARU,PARJ
COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
INTEGER MDCY,MDME,KFDP
DOUBLE PRECISION BRAT
COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
INTEGER MSEL,MSELPD,MSUB,KFIN
DOUBLE PRECISION CKIN
COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
INTEGER MSTP,MSTI
DOUBLE PRECISION PARP,PARI
COMMON/PYDATR/MRPY(6),RRPY(100)
INTEGER MRPY
DOUBLE PRECISION RRPY
C--use nuclear pdf?
COMMON/NPDF/MASS,NSET,EPS09,INITSTR
INTEGER NSET
DOUBLE PRECISION MASS
LOGICAL EPS09
CHARACTER*10 INITSTR
C--pdfset
common/pdf/pdfset
integer pdfset
C--number of protons
common/np/nproton
integer nproton
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--splitting integral
COMMON/SPLITINT/SPLITIGGV(1000,1000),SPLITIQQV(1000,1000),
&SPLITIQGV(1000,1000),QVAL(1000),ZMVAL(1000),QMAX,ZMMIN,NPOINT
INTEGER NPOINT
DOUBLE PRECISION SPLITIGGV,SPLITIQQV,SPLITIQGV,
&QVAL,ZMVAL,QMAX,ZMMIN
C--pdf common block
COMMON/PDFS/QINQX(2,1000),GINQX(2,1000),QINGX(2,1000),
&GINGX(2,1000)
DOUBLE PRECISION QINQX,GINQX,QINGX,GINGX
C--cross secttion common block
COMMON/XSECS/INTQ1(1001,101),INTQ2(1001,101),
&INTG1(1001,101),INTG2(1001,101)
DOUBLE PRECISION INTQ1,INTQ2,INTG1,INTG2
C--Sudakov common block
COMMON/INSUDA/SUDAQQ(1000,2),SUDAQG(1000,2),SUDAGG(1000,2)
&,SUDAGC(1000,2)
DOUBLE PRECISION SUDAQQ,SUDAQG,SUDAGG,SUDAGC
C--exponential integral for negative arguments
COMMON/EXPINT/EIX(3,1000),VALMAX,NVAL
INTEGER NVAL
DOUBLE PRECISION EIX,VALMAX
C--discard event flag
COMMON/DISC/NDISC,NSTRANGE,NGOOD,errcount,wdisc,DISCARD
LOGICAL DISCARD
INTEGER NDISC,NSTRANGE,NGOOD,errcount
double precision wdisc
C--factor in front of formation times
COMMON/FTIMEFAC/FTFAC
DOUBLE PRECISION FTFAC
C--factor in front of alphas argument
COMMON/ALPHASFAC/PTFAC
DOUBLE PRECISION PTFAC
C--number of scattering events
COMMON/CHECK/NSCAT,NSCATEFF,NSPLIT,nspliti,nsplitf,nistry,
&nisfail,nfsfail,nfstry,nttot,ntrej
DOUBLE PRECISION NSCAT,NSCATEFF,NSPLIT,nspliti,nsplitf,nistry,
&nisfail,nfsfail,nfstry,nttot,ntrej
C--number of extrapolations in tables
common/extrapolations/ntotspliti,noverspliti,ntotpdf,noverpdf,
&ntotxsec,noverxsec,ntotsuda,noversuda
integer ntotspliti,noverspliti,ntotpdf,noverpdf,
&ntotxsec,noverxsec,ntotsuda,noversuda
C--event weight
COMMON/WEIGHT/EVWEIGHT,sumofweights
double precision EVWEIGHT,sumofweights
C--event weight exponent
COMMON/WEXPO/WEIGHTEX
DOUBLE PRECISION WEIGHTEX
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
C--max rapidity
common/rapmax/etamax
double precision etamax
C--memory for error message from getdeltat
common/errline/errl
integer errl
C--organisation of event record
common/evrecord/nsim,npart,offset,hadrotype,sqrts,collider,hadro,
&shorthepmc,channel,isochannel
integer nsim,npart,offset,hadrotype
double precision sqrts
character*4 collider,channel
character*2 isochannel
logical hadro,shorthepmc
C--Pythia parameters
common/pythiaparams/PTMIN,PTMAX,weighted
double precision PTMIN,PTMAX
LOGICAL WEIGHTED
C--Variables local to this program
INTEGER NJOB,ios,pos,i,j,jj,intmass
DOUBLE PRECISION GETLTIMEMAX,EOVEST,r,pyr
character firstchar
CHARACTER*2 SNSET
CHARACTER*80 PDFFILE,XSECFILE,FILEMED,FILESPLIT,buffer,
&label,value
CHARACTER*100 HEPMCFILE,LOGFILE,FILENAME2
CHARACTER(LEN=100) filename
LOGICAL PDFEXIST,SPLITIEXIST,XSECEXIST
HPMCFID = 4
logfid = 3
C--default settings
nsim = 10000
njob = 0
logfile = 'out.log'
hepmcfile = 'out.hepmc'
filesplit = 'splitint.dat'
pdffile = 'pdfs.dat'
xsecfile = 'xsecs.dat'
filemed = 'medium-params.dat'
nf = 3
lqcd = 0.4
q0 = 1.5
ptmin = 5.
ptmax = 350.
etamax = 3.1
collider = 'PPJJ'
isochannel = 'XX'
channel = 'MUON'
sqrts = 2760
pdfset = 10042
nset = 1
mass = 208.
nproton = 82
weighted = .true.
weightex = 5.
angord = .true.
allhad = .false.
hadro = .true.
hadrotype = 0
shorthepmc = .true.
compress = .true.
lps = lqcd
scatrecoil = .false.
if (.not.hadro) shorthepmc = .true.
SCALEFACM=1.
ptfac=1.
ftfac=1.d0
if (iargc().eq.0) then
write(*,*)'No parameter file given, '//
&'will run with default settings.'
else
call getarg(1,filename)
write(*,*)'Reading parameters from ',filename
open(unit=1,file=filename,status='old',err=110)
do 120 i=1,1000
read(1, '(A)', iostat=ios) buffer
if(ios.ne.0) goto 130
firstchar = buffer(1:1)
if (firstchar.eq.'#') goto 120
pos=scan(buffer,' ')
label=buffer(1:pos)
value=buffer(pos+1:)
if(label.eq."NEVENT")then
read(value,*,iostat=ios) nsim
elseif(label.eq."NJOB")then
read(value,*,iostat=ios) njob
elseif(label.eq."LOGFILE")then
read(value,'(a)',iostat=ios) logfile
elseif(label.eq."HEPMCFILE")then
read(value,'(a)',iostat=ios) hepmcfile
elseif(label.eq."SPLITINTFILE")then
read(value,'(a)',iostat=ios) filesplit
elseif(label.eq."PDFFILE")then
read(value,'(a)',iostat=ios) pdffile
elseif(label.eq."XSECFILE")then
read(value,'(a)',iostat=ios) xsecfile
elseif(label.eq."MEDIUMPARAMS")then
read(value,'(a)',iostat=ios) filemed
elseif(label.eq."NF")then
read(value,*,iostat=ios) nf
elseif(label.eq."LAMBDAQCD")then
read(value,*,iostat=ios) lqcd
elseif(label.eq."Q0")then
read(value,*,iostat=ios) q0
elseif(label.eq."PTMIN")then
read(value,*,iostat=ios) ptmin
elseif(label.eq."PTMAX")then
read(value,*,iostat=ios) ptmax
elseif(label.eq."ETAMAX")then
read(value,*,iostat=ios) etamax
elseif(label.eq."PROCESS")then
read(value,*,iostat=ios) collider
elseif(label.eq."ISOCHANNEL")then
read(value,*,iostat=ios) isochannel
elseif(label.eq."CHANNEL")then
read(value,*,iostat=ios) channel
elseif(label.eq."SQRTS")then
read(value,*,iostat=ios) sqrts
elseif(label.eq."PDFSET")then
read(value,*,iostat=ios) pdfset
elseif(label.eq."NSET")then
read(value,*,iostat=ios) nset
elseif(label.eq."MASS")then
read(value,*,iostat=ios) mass
elseif(label.eq."NPROTON")then
read(value,*,iostat=ios) nproton
elseif(label.eq."WEIGHTED")then
read(value,*,iostat=ios) weighted
elseif(label.eq."WEXPO")then
read(value,*,iostat=ios) weightex
elseif(label.eq."ANGORD")then
read(value,*,iostat=ios) angord
elseif(label.eq."KEEPRECOILS")then
read(value,*,iostat=ios) allhad
elseif(label.eq."HADRO")then
read(value,*,iostat=ios) hadro
elseif(label.eq."HADROTYPE")then
read(value,*,iostat=ios) hadrotype
elseif(label.eq."SHORTHEPMC")then
read(value,*,iostat=ios) shorthepmc
elseif(label.eq."COMPRESS")then
read(value,*,iostat=ios) compress
else
write(*,*)'unknown label ',label
endif
120 continue
110 write(*,*)
& 'Unable to open parameter file, will exit the run.'
call exit(1)
130 close(1,status='keep')
write(*,*)'...done'
endif
if (ptmin.lt.3.d0) ptmin = 3.d0
OPEN(unit=logfid,file=LOGFILE,status='unknown')
MSTU(11)=logfid
call printtime
call printlogo(logfid)
write(logfid,*)
write(logfid,*)'parameters of the run:'
write(logfid,*)'NEVENT = ',nsim
write(logfid,*)'NJOB = ',njob
write(logfid,*)'LOGFILE = ',logfile
write(logfid,*)'HEPMCFILE = ',hepmcfile
write(logfid,*)'SPLITINTFILE = ',filesplit
write(logfid,*)'PDFFILE = ',pdffile
write(logfid,*)'XSECFILE = ',xsecfile
write(logfid,*)'MEDIUMPARAMS = ',filemed
write(logfid,*)'NF = ',nf
write(logfid,*)'LAMBDAQCD = ',lqcd
write(logfid,*)'Q0 = ',q0
write(logfid,*)'PTMIN = ',ptmin
write(logfid,*)'PTMAX = ',ptmax
write(logfid,*)'ETAMAX = ',etamax
write(logfid,*)'PROCESS = ',collider
write(logfid,*)'ISOCHANNEL = ',isochannel
write(logfid,*)'CHANNEL = ',channel
write(logfid,*)'SQRTS = ',sqrts
write(logfid,*)'PDFSET = ',pdfset
write(logfid,*)'NSET = ',nset
write(logfid,*)'MASS = ',mass
write(logfid,*)'NPROTON = ',nproton
write(logfid,*)'WEIGHTED = ',weighted
write(logfid,*)'WEXPO = ',weightex
write(logfid,*)'ANGORD = ',angord
write(logfid,*)'KEEPRECOILS = ',allhad
write(logfid,*)'HADRO = ',hadro
write(logfid,*)'HADROTYPE = ',hadrotype
write(logfid,*)'SHORTHEPMC = ',shorthepmc
write(logfid,*)'COMPRESS = ',compress
write(logfid,*)
call flush(logfid)
if ((collider.ne.'PPJJ').and.(collider.ne.'EEJJ')
& .and.(collider.ne.'PPYJ').and.(collider.ne.'PPYQ')
& .and.(collider.ne.'PPYG')
& .and.(collider.ne.'PPZJ').and.(collider.ne.'PPZQ')
& .and.(collider.ne.'PPZG').and.(collider.ne.'PPWJ')
& .and.(collider.ne.'PPWQ').and.(collider.ne.'PPWG')
& .and.(collider.ne.'PPDY')) then
write(logfid,*)'Fatal error: colliding system unknown, '//
& 'will exit now'
call exit(1)
endif
C--initialize medium
intmass = int(mass)
CALL MEDINIT(FILEMED,logfid,etamax,intmass)
CALL MEDNEXTEVT
OPEN(unit=HPMCFID,file=HEPMCFILE,status='unknown')
WRITE(HPMCFID,*)
WRITE(HPMCFID,'(A)')'HepMC::Version 2.06.05'
WRITE(HPMCFID,'(A)')'HepMC::IO_GenEvent-START_EVENT_LISTING'
NPART=2
if(ptmax.gt.0.)then
EOVEST=MIN(1.5*(PTMAX+50.)*COSH(ETAMAX),sqrts/2.)
else
EOVEST=sqrts/2.
endif
CALL EIXINT
CALL INSUDAINT(EOVEST)
write(logfid,*)
INQUIRE(file=FILESPLIT,exist=SPLITIEXIST)
IF(SPLITIEXIST)THEN
write(logfid,*)'read splitting integrals from ',FILESPLIT
OPEN(unit=10,file=FILESPLIT,status='old')
READ(10,*)QMAX,ZMMIN,NPOINT
DO 893 I=1,NPOINT+1
READ(10,*) QVAL(I),ZMVAL(I)
893 CONTINUE
DO 891 I=1,NPOINT+1
DO 892 J=1,NPOINT+1
READ(10,*)SPLITIGGV(I,J),SPLITIQQV(I,J),SPLITIQGV(I,J)
892 CONTINUE
891 CONTINUE
CLOSE(10,status='keep')
ELSE
write(logfid,*)'have to integrate splitting functions, '//
&'this may take some time'
CALL SPLITFNCINT(EOVEST)
INQUIRE(file=FILESPLIT,exist=SPLITIEXIST)
IF(.NOT.SPLITIEXIST)THEN
write(logfid,*)'write splitting integrals to ',FILESPLIT
OPEN(unit=10,file=FILESPLIT,status='new')
WRITE(10,*)QMAX,ZMMIN,NPOINT
DO 896 I=1,NPOINT+1
WRITE(10,*) QVAL(I),ZMVAL(I)
896 CONTINUE
DO 897 I=1,NPOINT+1
DO 898 J=1,NPOINT+1
WRITE(10,*)SPLITIGGV(I,J),SPLITIQQV(I,J),SPLITIQGV(I,J)
898 CONTINUE
897 CONTINUE
CLOSE(10,status='keep')
ENDIF
ENDIF
write(logfid,*)
INQUIRE(file=PDFFILE,exist=PDFEXIST)
IF(PDFEXIST)THEN
write(logfid,*)'read pdfs from ',PDFFILE
OPEN(unit=10,file=PDFFILE,status='old')
DO 872 I=1,2
DO 873 J=1,1000
READ(10,*)QINQX(I,J),GINQX(I,J),QINGX(I,J),GINGX(I,J)
873 CONTINUE
872 CONTINUE
CLOSE(10,status='keep')
ELSE
write(logfid,*)'have to integrate pdfs, this may take some time'
CALL PDFINT(EOVEST)
INQUIRE(file=PDFFILE,exist=PDFEXIST)
IF(.NOT.PDFEXIST)THEN
write(logfid,*)'write pdfs to ',PDFFILE
OPEN(unit=10,file=PDFFILE,status='new')
DO 876 I=1,2
DO 877 J=1,1000
WRITE(10,*)QINQX(I,J),GINQX(I,J),QINGX(I,J),GINGX(I,J)
877 CONTINUE
876 CONTINUE
CLOSE(10,status='keep')
ENDIF
ENDIF
write(logfid,*)
INQUIRE(file=XSECFILE,exist=XSECEXIST)
IF(XSECEXIST)THEN
write(logfid,*)'read cross sections from ',XSECFILE
OPEN(unit=10,file=XSECFILE,status='old')
DO 881 J=1,1001
DO 885 JJ=1,101
READ(10,*)INTQ1(J,JJ),INTQ2(J,JJ),
&INTG1(J,JJ),INTG2(J,JJ)
885 CONTINUE
881 CONTINUE
CLOSE(10,status='keep')
ELSE
write(logfid,*)'have to integrate cross sections, '//
&'this may take some time'
CALL XSECINT(EOVEST)
INQUIRE(file=XSECFILE,exist=XSECEXIST)
IF(.NOT.XSECEXIST)THEN
write(logfid,*)'write cross sections to ',XSECFILE
OPEN(unit=10,file=XSECFILE,status='new')
DO 883 J=1,1001
DO 884 JJ=1,101
WRITE(10,*)INTQ1(J,JJ),INTQ2(J,JJ),
&INTG1(J,JJ),INTG2(J,JJ)
884 CONTINUE
883 CONTINUE
CLOSE(10,status='keep')
ENDIF
ENDIF
write(logfid,*)
CALL FLUSH(3)
C--initialise random number generator status
IF(NJOB.GT.0)THEN
MRPY(1)=NJOB*1000
MRPY(2)=0
ENDIF
C--Call PYR once for initialization
R=PYR(0)
NDISC=0
NGOOD=0
NSTRANGE=0
ERRCOUNT=0
errl = 0
NSCAT=0.d0
NSCATEFF=0.d0
NSPLIT=0.d0
nspliti=0.d0
nsplitf=0.d0
nistry=0.d0
nisfail=0.d0
nfstry=0.d0
nfsfail=0.d0
nttot=0.d0
ntrej=0.d0
ntotspliti=0
noverspliti=0
ntotpdf=0
noverpdf=0
ntotxsec=0
noverxsec=0
ntotsuda=0
noversuda=0
IF(NSET.EQ.0)THEN
EPS09=.FALSE.
ELSE
EPS09=.TRUE.
IF(NSET.LT.10)THEN
WRITE(SNSET,'(i1)') NSET
ELSE
WRITE(SNSET,'(i2)') NSET
ENDIF
INITSTR='EPS09LO,'//SNSET
ENDIF
end
***********************************************************************
*** subroutine initpythia
***********************************************************************
subroutine initpythia(beam1,beam2)
implicit none
INTEGER PYCOMP
INTEGER NMXHEP
C--Common block of Pythia
COMMON/PYJETS/N,NPAD,K(23000,5),P(23000,5),V(23000,5)
INTEGER N,NPAD,K
DOUBLE PRECISION P,V
COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
INTEGER MSTU,MSTJ
DOUBLE PRECISION PARU,PARJ
COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
INTEGER MDCY,MDME,KFDP
DOUBLE PRECISION BRAT
COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
INTEGER MSEL,MSELPD,MSUB,KFIN
DOUBLE PRECISION CKIN
COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
INTEGER MSTP,MSTI
DOUBLE PRECISION PARP,PARI
COMMON/PYDATR/MRPY(6),RRPY(100)
INTEGER MRPY
DOUBLE PRECISION RRPY
C--use nuclear pdf?
COMMON/NPDF/MASS,NSET,EPS09,INITSTR
INTEGER NSET
DOUBLE PRECISION MASS
LOGICAL EPS09
CHARACTER*10 INITSTR
C--pdfset
common/pdf/pdfset
integer pdfset
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--discard event flag
COMMON/DISC/NDISC,NSTRANGE,NGOOD,errcount,wdisc,DISCARD
LOGICAL DISCARD
INTEGER NDISC,NSTRANGE,NGOOD,errcount
double precision wdisc
C--event weight
COMMON/WEIGHT/EVWEIGHT,sumofweights
double precision EVWEIGHT,sumofweights
C--event weight exponent
COMMON/WEXPO/WEIGHTEX
DOUBLE PRECISION WEIGHTEX
C--memory for error message from getdeltat
common/errline/errl
integer errl
C--organisation of event record
common/evrecord/nsim,npart,offset,hadrotype,sqrts,collider,hadro,
&shorthepmc,channel,isochannel
integer nsim,npart,offset,hadrotype
double precision sqrts
character*4 collider,channel
character*2 isochannel
logical hadro,shorthepmc
C--Pythia parameters
common/pythiaparams/PTMIN,PTMAX,weighted
double precision PTMIN,PTMAX
LOGICAL WEIGHTED
C--Variables local to this program
character*2 beam1,beam2
C--initialise PYTHIA
C--no multiple interactions
MSTP(81) = 0
C--initial state radiation
MSTP(61)=1
C--switch off final state radiation
MSTP(71)=0
C--No hadronisation (yet)
MSTP(111)=0
C--parameter affecting treatment of string corners
PARU(14)=1.
C--Min shat in simulation
CKIN(1)=2.
C--pT-cut
CKIN(3)=PTMIN
CKIN(4)=PTMAX
C--use LHAPDF
MSTP(52)=2
C--choose pdf: CTEQ6ll (LO fit/LO alphas) - 10042
C MSTW2008 (LO central) - 21000
MSTP(51)=PDFSET
IF(COLLIDER.EQ.'PPYQ')THEN
MSEL=0
MSUB(29)=1
ELSEIF(COLLIDER.EQ.'PPYG')THEN
MSEL=0
MSUB(14)=1
MSUB(115)=1
ELSEIF(COLLIDER.EQ.'PPYJ')THEN
MSEL=0
MSUB(14)=1
MSUB(29)=1
MSUB(115)=1
ELSEIF((COLLIDER.EQ.'PPZJ').or.(COLLIDER.EQ.'PPZQ')
& .or.(COLLIDER.EQ.'PPZG')
& .or.(collider.eq.'PPDY'))THEN
MSEL=0
IF((COLLIDER.EQ.'PPZJ').or.(COLLIDER.EQ.'PPZQ')) MSUB(30)=1
IF((COLLIDER.EQ.'PPZJ').or.(COLLIDER.EQ.'PPZG')) MSUB(15)=1
IF(COLLIDER.EQ.'PPDY') MSUB(1)=1
MDME(174,1)=0 !Z decay into d dbar',
MDME(175,1)=0 !Z decay into u ubar',
MDME(176,1)=0 !Z decay into s sbar',
MDME(177,1)=0 !Z decay into c cbar',
MDME(178,1)=0 !Z decay into b bbar',
MDME(179,1)=0 !Z decay into t tbar',
MDME(182,1)=0 !Z decay into e- e+',
MDME(183,1)=0 !Z decay into nu_e nu_ebar',
MDME(184,1)=0 !Z decay into mu- mu+',
MDME(185,1)=0 !Z decay into nu_mu nu_mubar',
MDME(186,1)=0 !Z decay into tau- tau+',
MDME(187,1)=0 !Z decay into nu_tau nu_taubar',
if (channel.EQ.'ELEC')THEN
MDME(182,1)=1
ELSEIF(channel.EQ.'MUON')THEN
MDME(184,1)=1
ENDIF
ELSEIF((COLLIDER.EQ.'PPWJ').or.(COLLIDER.EQ.'PPWQ')
& .or.(COLLIDER.EQ.'PPWG'))THEN
MSEL=0
IF((COLLIDER.EQ.'PPWJ').or.(COLLIDER.EQ.'PPWQ')) MSUB(31)=1
IF((COLLIDER.EQ.'PPWJ').or.(COLLIDER.EQ.'PPWG')) MSUB(16)=1
MDME(190,1)=0 ! W+ decay into dbar u,
MDME(191,1)=0 ! W+ decay into dbar c,
MDME(192,1)=0 ! W+ decay into dbar t,
MDME(194,1)=0 ! W+ decay into sbar u,
MDME(195,1)=0 ! W+ decay into sbar c,
MDME(196,1)=0 ! W+ decay into sbar t,
MDME(198,1)=0 ! W+ decay into bbar u,
MDME(199,1)=0 ! W+ decay into bbar c,
MDME(200,1)=0 ! W+ decay into bbar t,
MDME(202,1)=0 ! W+ decay into b'bar u,
MDME(203,1)=0 ! W+ decay into b'bar c,
MDME(204,1)=0 ! W+ decay into b'bar t,
MDME(206,1)=0 ! W+ decay into e+ nu_e,
MDME(207,1)=0 ! W+ decay into mu+ nu_mu,
MDME(208,1)=0 ! W+ decay into tau+ nu_tau,
MDME(209,1)=0 ! W+ decay into tau'+ nu'_tau,
if (channel.EQ.'ELEC')THEN
MDME(206,1)=1
ELSEIF(channel.EQ.'MUON')THEN
MDME(207,1)=1
ENDIF
ELSE
C--All QCD processes are active
MSEL=1
ENDIF
! MSEL=0
! MSUB(11)=1
! MSUB(12)=1
! MSUB(53)=1
! MSUB(13)=1
! MSUB(68)=1
! MSUB(28)=1
C--weighted events
IF(WEIGHTED) MSTP(142)=1
C--number of errors to be printed
MSTU(22)=MAX(10,INT(5.*NSIM/100.))
C--number of lines in event record
MSTU(4)=23000
MSTU(5)=23000
C--switch off pi0 decay
MDCY(PYCOMP(111),1)=0
C--initialisation call
IF(COLLIDER.EQ.'EEJJ')THEN
OFFSET=9
CALL PYINIT('CMS',beam1,beam2,sqrts)
ELSEIF((COLLIDER.EQ.'PPJJ').OR.(COLLIDER.EQ.'PPYJ').OR.
& (COLLIDER.EQ.'PPYG').OR.(COLLIDER.EQ.'PPYQ'))THEN
OFFSET=8
CALL PYINIT('CMS',beam1,beam2,sqrts)
ELSEIF((COLLIDER.EQ.'PPWJ').OR.(COLLIDER.EQ.'PPZJ').or.
& (COLLIDER.EQ.'PPWQ').OR.(COLLIDER.EQ.'PPZQ').or.
& (COLLIDER.EQ.'PPWG').OR.(COLLIDER.EQ.'PPZG'))THEN
OFFSET=10
CALL PYINIT('CMS',beam1,beam2,sqrts)
elseif (collider.eq.'PPDY') then
CALL PYINIT('CMS',beam1,beam2,sqrts)
ENDIF
end
***********************************************************************
*** subroutine genevent
***********************************************************************
subroutine genevent(j,b1,b2)
implicit none
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
INTEGER PYCOMP
INTEGER NMXHEP
C--Common block of Pythia
COMMON/PYJETS/N,NPAD,K(23000,5),P(23000,5),V(23000,5)
INTEGER N,NPAD,K
DOUBLE PRECISION P,V
COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
INTEGER MSTU,MSTJ
DOUBLE PRECISION PARU,PARJ
COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
INTEGER MDCY,MDME,KFDP
DOUBLE PRECISION BRAT
COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
INTEGER MSEL,MSELPD,MSUB,KFIN
DOUBLE PRECISION CKIN
COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
INTEGER MSTP,MSTI
DOUBLE PRECISION PARP,PARI
COMMON/PYDATR/MRPY(6),RRPY(100)
INTEGER MRPY
DOUBLE PRECISION RRPY
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--discard event flag
COMMON/DISC/NDISC,NSTRANGE,NGOOD,errcount,wdisc,DISCARD
LOGICAL DISCARD
INTEGER NDISC,NSTRANGE,NGOOD,errcount
double precision wdisc
C--variables for angular ordering
COMMON/ANGOR/ZA(23000),ZD(23000),THETAA(23000),QQBARD(23000)
DOUBLE PRECISION ZA,ZD,THETAA
LOGICAL QQBARD
C--factor in front of formation times
COMMON/FTIMEFAC/FTFAC
DOUBLE PRECISION FTFAC
C--time common block
COMMON/TIME/MV(23000,5)
DOUBLE PRECISION MV
C--colour index common block
COMMON/COLOUR/TRIP(23000),ANTI(23000),COLMAX
INTEGER TRIP,ANTI,COLMAX
C--number of scattering events
COMMON/CHECK/NSCAT,NSCATEFF,NSPLIT,nspliti,nsplitf,nistry,
&nisfail,nfsfail,nfstry,nttot,ntrej
DOUBLE PRECISION NSCAT,NSCATEFF,NSPLIT,nspliti,nsplitf,nistry,
&nisfail,nfsfail,nfstry,nttot,ntrej
C--event weight
COMMON/WEIGHT/EVWEIGHT,sumofweights
double precision EVWEIGHT,sumofweights
C--event weight exponent
COMMON/WEXPO/WEIGHTEX
DOUBLE PRECISION WEIGHTEX
C--max rapidity
common/rapmax/etamax
double precision etamax
C--production point
common/jetpoint/x0,y0
double precision x0,y0
C--initial pt and virtuality
common/injet/isgluon(2),inpt(2),inmass(2),inphi(2),ineta(2)
integer isgluon
double precision inpt,inmass,inphi,ineta
C--organisation of event record
common/evrecord/nsim,npart,offset,hadrotype,sqrts,collider,hadro,
&shorthepmc,channel,isochannel
integer nsim,npart,offset,hadrotype
double precision sqrts
character*4 collider,channel
character*2 isochannel
logical hadro,shorthepmc
C--Variables local to this program
INTEGER NOLD,PID,IPART,LME1,LME2,j,i,LME1ORIG,LME2ORIG,llep1,
&llep2,lv
DOUBLE PRECISION PYR,ENI,QMAX1,R,GETMASS,PYP,Q1,Q2,P21,P22,ETOT,
&QMAX2,POLD,EN1,EN2,BETA(3),ENEW1,ENEW2,emax,lambda,x1,x2,x3,
&MEWEIGHT,PSWEIGHT,WEIGHT,EPS1,EPS2,THETA1,THETA2,Z1,Z2,
&getltimemax,pi,m1,m2
character*2 b1,b2
LOGICAL FIRSTTRIP,WHICH1,WHICH2,ISDIQUARK
DATA PI/3.141592653589793d0/
N=0
COLMAX=600
DISCARD=.FALSE.
DO 91 I=1,23000
MV(I,1)=0.d0
MV(I,2)=0.d0
MV(I,3)=0.d0
MV(I,4)=0.d0
MV(I,5)=0.d0
91 CONTINUE
CALL MEDNEXTEVT
C--initialisation with matrix element
C--production vertex
CALL PICKVTX(X0,Y0)
LTIME=GETLTIMEMAX()
99 CALL PYEVNT
NPART=N-OFFSET
EVWEIGHT=PARI(10)
SUMOFWEIGHTS=SUMOFWEIGHTS+EVWEIGHT
IF((COLLIDER.EQ.'EEJJ').AND.(ABS(K(8,2)).GT.6))THEN
WDISC=WDISC+EVWEIGHT
NDISC=NDISC+1
GOTO 102
ELSE
NGOOD=NGOOD+1
ENDIF
C--DY: don't have to do anything
if (collider.eq.'PPDY') then
CALL PYEXEC
call CONVERTTOHEPMC(HPMCFID,NGOOD,PID,b1,b2)
goto 102
endif
C-- prepare event record
if((COLLIDER.EQ.'PPZJ').OR.(COLLIDER.EQ.'PPZQ').or.
& (COLLIDER.EQ.'PPZG').or.(COLLIDER.EQ.'PPWJ').or.
& (COLLIDER.EQ.'PPWQ').or.(COLLIDER.EQ.'PPWG'))THEN
LME1ORIG=7
LME2ORIG=8
if(abs(k(7,2)).gt.21) then
lv=7
else
lv=8
endif
ELSE
LME1ORIG=OFFSET-1
LME2ORIG=OFFSET
ENDIF
DO 180 IPART=OFFSET+1, OFFSET+NPART
C--find decay leptons in V+jet events
if((COLLIDER.EQ.'PPZJ').OR.(COLLIDER.EQ.'PPZQ').or.
& (COLLIDER.EQ.'PPZG').or.(COLLIDER.EQ.'PPWJ').or.
& (COLLIDER.EQ.'PPWQ').or.(COLLIDER.EQ.'PPWG'))THEN
if(k(ipart,3).eq.offset-1) llep1=ipart
if(k(ipart,3).eq.offset) llep2=ipart
endif
IF(K(IPART,3).EQ.(LME1ORIG))THEN
LME1=IPART
ELSEIF(K(IPART,3).EQ.LME2ORIG)THEN
LME2=IPART
ELSE
TRIP(IPART)=0
ANTI(IPART)=0
ZD(IPART)=0.d0
THETAA(IPART)=0.d0
ENDIF
C--assign colour indices
IF(K(IPART,1).EQ.2)THEN
IF(K(IPART-1,1).EQ.2)THEN
C--in middle of colour singlet
IF(FIRSTTRIP)THEN
TRIP(IPART)=COLMAX+1
ANTI(IPART)=TRIP(IPART-1)
ELSE
TRIP(IPART)=ANTI(IPART-1)
ANTI(IPART)=COLMAX+1
ENDIF
COLMAX=COLMAX+1
ELSE
C--beginning of colour singlet
IF(((ABS(K(IPART,2)).LT.10).AND.(K(IPART,2).GT.0))
& .OR.(ISDIQUARK(K(IPART,2)).AND.(K(IPART,2).LT.0)))THEN
TRIP(IPART)=COLMAX+1
ANTI(IPART)=0
FIRSTTRIP=.TRUE.
ELSE
TRIP(IPART)=0
ANTI(IPART)=COLMAX+1
FIRSTTRIP=.FALSE.
ENDIF
COLMAX=COLMAX+1
ENDIF
ENDIF
IF(K(IPART,1).EQ.1)THEN
C--end of colour singlet
IF(FIRSTTRIP)THEN
TRIP(IPART)=0
ANTI(IPART)=TRIP(IPART-1)
ELSE
TRIP(IPART)=ANTI(IPART-1)
ANTI(IPART)=0
ENDIF
ENDIF
180 CONTINUE
if (k(lme1,1).lt.11) K(LME1,1)=1
if (k(lme2,1).lt.11) K(LME2,1)=1
PID=K(LME1,2)
ENI=MAX(P(LME1,4),P(LME2,4))
DO 183 IPART=OFFSET+1, OFFSET+NPART
IF((IPART.NE.LME1).AND.(IPART.NE.LME2).AND.(K(IPART,1).LT.11))
& K(IPART,1)=4
if (k(ipart,2).eq.22) k(ipart,1)=4
183 CONTINUE
C--find virtualities and adapt four-vectors
if((COLLIDER.EQ.'PPZJ').OR.(COLLIDER.EQ.'PPZQ').or.
& (COLLIDER.EQ.'PPZG').or.(COLLIDER.EQ.'PPWJ').or.
& (COLLIDER.EQ.'PPWQ').or.(COLLIDER.EQ.'PPWG'))THEN
if (abs(k(lme1,2)).gt.21) then
QMAX1=0.d0
QMAX2=sqrt(pari(18)+p(lme1,5)**2)
else
QMAX1=sqrt(pari(18)+p(lme2,5)**2)
QMAX2=0.d0
endif
EMAX=P(LME1,4)+P(LME2,4)
THETA1=-1.d0
THETA2=-1.d0
ELSEIF(COLLIDER.EQ.'PPJJ'.OR.COLLIDER.EQ.'PPYJ'
& .OR.COLLIDER.EQ.'PPYQ'.OR.COLLIDER.EQ.'PPYG')THEN
if (k(lme1,1).eq.4) then
qmax1 = 0.d0
else
QMAX1=pari(17)
endif
if (k(lme2,1).eq.4) then
qmax2 = 0.d0
else
QMAX2=pari(17)
endif
! QMAX1=PYP(LME1,10)*exp(0.3*abs(pyp(lme1,17)-pyp(lme2,17))/2.)/2.
! QMAX2=PYP(LME2,10)*exp(0.3*abs(pyp(lme1,17)-pyp(lme2,17))/2.)/2.
EMAX=P(LME1,4)+P(LME2,4)
THETA1=-1.d0
THETA2=-1.d0
ENDIF
EN1=P(LME1,4)
EN2=P(LME2,4)
BETA(1)=(P(LME1,1)+P(LME2,1))/(P(LME1,4)+P(LME2,4))
BETA(2)=(P(LME1,2)+P(LME2,2))/(P(LME1,4)+P(LME2,4))
BETA(3)=(P(LME1,3)+P(LME2,3))/(P(LME1,4)+P(LME2,4))
CALL PYROBO(LME1,LME1,0d0,0d0,-BETA(1),-BETA(2),-BETA(3))
CALL PYROBO(LME2,LME2,0d0,0d0,-BETA(1),-BETA(2),-BETA(3))
ETOT=P(LME1,4)+P(LME2,4)
IF(COLLIDER.EQ.'EEJJ')THEN
QMAX1=ETOT
QMAX2=ETOT
EMAX=P(LME1,4)+P(LME2,4)
THETA1=-1.d0
THETA2=-1.d0
ENDIF
C-- find virtuality
Q1=GETMASS(0.d0,QMAX1,THETA1,EMAX,k(lme1,2),EMAX,.FALSE.,
& Z1,WHICH1)
Q2=GETMASS(0.d0,QMAX2,THETA2,EMAX,k(lme2,2),EMAX,.FALSE.,
& Z2,WHICH2)
182 if (abs(k(lme1,2)).gt.21) then
m1=p(lme1,5)
else
m1=q1
endif
if (abs(k(lme2,2)).gt.21) then
m2=p(lme2,5)
else
m2=q2
endif
ENEW1=ETOT/2.d0 + (m1**2-m2**2)/(2.*ETOT)
ENEW2=ETOT/2.d0 - (m1**2-m2**2)/(2.*ETOT)
P21 = (ETOT/2.d0 + (m1**2-m2**2)/(2.*ETOT))**2 - m1**2
P22 = (ETOT/2.d0 - (m1**2-m2**2)/(2.*ETOT))**2 - m2**2
WEIGHT=1.d0
IF((PYR(0).GT.WEIGHT).OR.(P21.LT.0.d0).OR.(P22.LT.0.d0)
& .OR.(ENEW1.LT.0.d0).OR.(ENEW2.LT.0.d0)
& )THEN
IF(Q1.GT.Q2)THEN
Q1=GETMASS(0.d0,Q1,THETA1,EMAX,k(lme1,2),EMAX,.FALSE.,
& Z1,WHICH1)
ELSE
Q2=GETMASS(0.d0,Q2,THETA2,EMAX,k(lme2,2),EMAX,.FALSE.,
& Z2,WHICH2)
ENDIF
GOTO 182
ENDIF
POLD=PYP(LME1,8)
P(LME1,1)=P(LME1,1)*SQRT(P21)/POLD
P(LME1,2)=P(LME1,2)*SQRT(P21)/POLD
P(LME1,3)=P(LME1,3)*SQRT(P21)/POLD
P(LME1,4)=ENEW1
P(LME1,5)=m1
POLD=PYP(LME2,8)
P(LME2,1)=P(LME2,1)*SQRT(P22)/POLD
P(LME2,2)=P(LME2,2)*SQRT(P22)/POLD
P(LME2,3)=P(LME2,3)*SQRT(P22)/POLD
P(LME2,4)=ENEW2
P(LME2,5)=m2
CALL PYROBO(LME1,LME1,0d0,0d0,BETA(1),BETA(2),BETA(3))
CALL PYROBO(LME2,LME2,0d0,0d0,BETA(1),BETA(2),BETA(3))
C--correct for overestimated energy
IF(Q1.GT.0.d0)THEN
EPS1=0.5-0.5*SQRT(1.-Q0**2/Q1**2)
& *SQRT(1.-Q1**2/P(LME1,4)**2)
IF((Z1.LT.EPS1).OR.(Z1.GT.(1.-EPS1)))THEN
Q1=GETMASS(0.d0,Q1,THETA1,EMAX,k(lme1,2),EMAX,.FALSE.,
& Z1,WHICH1)
CALL PYROBO(LME1,LME1,0d0,0d0,-BETA(1),-BETA(2),-BETA(3))
CALL PYROBO(LME2,LME2,0d0,0d0,-BETA(1),-BETA(2),-BETA(3))
GOTO 182
ENDIF
ENDIF
IF(Q2.GT.0.d0)THEN
EPS2=0.5-0.5*SQRT(1.-Q0**2/Q2**2)
& *SQRT(1.-Q2**2/P(LME2,4)**2)
IF((Z2.LT.EPS2).OR.(Z2.GT.(1.-EPS2)))THEN
Q2=GETMASS(0.d0,Q2,THETA2,EMAX,k(lme2,2),EMAX,.FALSE.,
& Z2,WHICH2)
CALL PYROBO(LME1,LME1,0d0,0d0,-BETA(1),-BETA(2),-BETA(3))
CALL PYROBO(LME2,LME2,0d0,0d0,-BETA(1),-BETA(2),-BETA(3))
GOTO 182
ENDIF
ENDIF
C--correct to ME for first parton
IF(COLLIDER.EQ.'EEJJ')THEN
BETA(1)=(P(LME1,1)+P(LME2,1))/(P(LME1,4)+P(LME2,4))
BETA(2)=(P(LME1,2)+P(LME2,2))/(P(LME1,4)+P(LME2,4))
BETA(3)=(P(LME1,3)+P(LME2,3))/(P(LME1,4)+P(LME2,4))
CALL PYROBO(LME1,LME1,0d0,0d0,-BETA(1),-BETA(2),-BETA(3))
CALL PYROBO(LME2,LME2,0d0,0d0,-BETA(1),-BETA(2),-BETA(3))
IF(Q1.GT.0.d0)THEN
C--generate z value
X1=Z1*(ETOT**2+Q1**2)/ETOT**2
X2=(ETOT**2-Q1**2)/ETOT**2
X3=(1.-Z1)*(ETOT**2+Q1**2)/ETOT**2
PSWEIGHT=(1.-X1)*(1.+(X1/(2.-X2))**2)/X3
& + (1.-X2)*(1.+(X2/(2.-X1))**2)/X3
MEWEIGHT=X1**2+X2**2
WEIGHT=MEWEIGHT/PSWEIGHT
IF(PYR(0).GT.WEIGHT)THEN
184 Q1=GETMASS(0.d0,Q1,THETA1,EMAX,k(lme1,2),EMAX,.FALSE.,
& Z1,WHICH1)
ENDIF
ENDIF
C--correct to ME for second parton
IF(Q2.GT.0.d0)THEN
C--generate z value
X1=(ETOT**2-Q2**2)/ETOT**2
X2=Z2*(ETOT**2+Q2**2)/ETOT**2
X3=(1.-Z2)*(ETOT**2+Q2**2)/ETOT**2
PSWEIGHT=(1.-X1)*(1.+(X1/(2.-X2))**2)/X3
& + (1.-X2)*(1.+(X2/(2.-X1))**2)/X3
MEWEIGHT=X1**2+X2**2
WEIGHT=MEWEIGHT/PSWEIGHT
IF(PYR(0).GT.WEIGHT)THEN
185 Q2=GETMASS(0.d0,Q2,THETA2,EMAX,k(lme1,2),EMAX,.FALSE.,
& Z2,WHICH2)
ENDIF
ENDIF
186 ENEW1=ETOT/2.d0 + (Q1**2-Q2**2)/(2.*ETOT)
ENEW2=ETOT/2.d0 - (Q1**2-Q2**2)/(2.*ETOT)
P21 = (ETOT/2.d0 + (Q1**2-Q2**2)/(2.*ETOT))**2 - Q1**2
P22 = (ETOT/2.d0 - (Q1**2-Q2**2)/(2.*ETOT))**2 - Q2**2
POLD=PYP(LME1,8)
P(LME1,1)=P(LME1,1)*SQRT(P21)/POLD
P(LME1,2)=P(LME1,2)*SQRT(P21)/POLD
P(LME1,3)=P(LME1,3)*SQRT(P21)/POLD
P(LME1,4)=ENEW1
P(LME1,5)=Q1
POLD=PYP(LME2,8)
P(LME2,1)=P(LME2,1)*SQRT(P22)/POLD
P(LME2,2)=P(LME2,2)*SQRT(P22)/POLD
P(LME2,3)=P(LME2,3)*SQRT(P22)/POLD
P(LME2,4)=ENEW2
P(LME2,5)=Q2
CALL PYROBO(LME1,LME1,0d0,0d0,BETA(1),BETA(2),BETA(3))
CALL PYROBO(LME2,LME2,0d0,0d0,BETA(1),BETA(2),BETA(3))
C--correct for overestimated energy
IF(Q1.GT.0.d0)THEN
EPS1=0.5-0.5*SQRT(1.-Q0**2/Q1**2)
& *SQRT(1.-Q1**2/P(LME1,4)**2)
IF((Z1.LT.EPS1).OR.(Z1.GT.(1.-EPS1)))THEN
Q1=GETMASS(0.d0,Q1,THETA1,EMAX,k(lme1,2),EMAX,.FALSE.,
& Z1,WHICH1)
CALL PYROBO(LME1,LME1,0d0,0d0,-BETA(1),-BETA(2),-BETA(3))
CALL PYROBO(LME2,LME2,0d0,0d0,-BETA(1),-BETA(2),-BETA(3))
GOTO 186
ENDIF
ENDIF
IF(Q2.GT.0.d0)THEN
EPS2=0.5-0.5*SQRT(1.-Q0**2/Q2**2)
& *SQRT(1.-Q2**2/P(LME2,4)**2)
IF((Z2.LT.EPS2).OR.(Z2.GT.(1.-EPS2)))THEN
Q2=GETMASS(0.d0,Q2,THETA2,EMAX,k(lme2,2),EMAX,.FALSE.,
& Z2,WHICH2)
CALL PYROBO(LME1,LME1,0d0,0d0,-BETA(1),-BETA(2),-BETA(3))
CALL PYROBO(LME2,LME2,0d0,0d0,-BETA(1),-BETA(2),-BETA(3))
GOTO 186
ENDIF
ENDIF
ENDIF
C--transfer recoil to decay leptons in V+jet
if((COLLIDER.EQ.'PPZJ').OR.(COLLIDER.EQ.'PPZQ').or.
& (COLLIDER.EQ.'PPZG').or.(COLLIDER.EQ.'PPWJ').or.
& (COLLIDER.EQ.'PPWQ').or.(COLLIDER.EQ.'PPWG'))THEN
beta(1)=p(lv,1)/p(lv,4)
beta(2)=p(lv,2)/p(lv,4)
beta(3)=p(lv,3)/p(lv,4)
CALL PYROBO(llep1,llep1,0d0,0d0,-BETA(1),-BETA(2),-BETA(3))
CALL PYROBO(llep2,llep2,0d0,0d0,-BETA(1),-BETA(2),-BETA(3))
if (abs(k(lme1,2)).gt.21) then
beta(1)=p(lme1,1)/p(lme1,4)
beta(2)=p(lme1,2)/p(lme1,4)
beta(3)=p(lme1,3)/p(lme1,4)
else
beta(1)=p(lme2,1)/p(lme2,4)
beta(2)=p(lme2,2)/p(lme2,4)
beta(3)=p(lme2,3)/p(lme2,4)
endif
CALL PYROBO(llep1,llep1,0d0,0d0,BETA(1),BETA(2),BETA(3))
CALL PYROBO(llep2,llep2,0d0,0d0,BETA(1),BETA(2),BETA(3))
endif
C--store initial parton pt and mass for output
if (k(lme1,1).eq.1) then
inpt(1) = pyp(lme1,10)
! inpt(1) = p(lme1,4)*sin(pyp(lme1,13))
inmass(1) = p(lme1,5)
inphi(1) = pyp(lme1,15)
ineta(1) = pyp(lme1,19)
inpt(2) = pyp(lme2,10)
! inpt(2) = p(lme2,4)*sin(pyp(lme2,13))
inmass(2) = p(lme2,5)
inphi(2) = pyp(lme2,15)
ineta(2) = pyp(lme2,19)
if (k(lme1,2).eq.21) then
isgluon(1) = 1
elseif (abs(k(lme1,2)).le.5) then
isgluon(1) = 0
else
isgluon(1) = 2
endif
if (k(lme2,2).eq.21) then
isgluon(2) = 1
elseif (abs(k(lme2,2)).le.5) then
isgluon(2) = 0
else
isgluon(2) = 2
endif
else
inpt(1) = pyp(lme2,10)
! inpt(1) = p(lme2,4)*sin(pyp(lme2,13))
inmass(1) = p(lme2,5)
inphi(1) = pyp(lme2,15)
ineta(1) = pyp(lme2,19)
inpt(2) = pyp(lme1,10)
! inpt(2) = p(lme1,4)*sin(pyp(lme1,13))
inmass(2) = p(lme1,5)
inphi(2) = pyp(lme1,15)
ineta(2) = pyp(lme1,19)
if (k(lme2,2).eq.21) then
isgluon(1) = 1
elseif (abs(k(lme2,2)).le.5) then
isgluon(1) = 0
else
isgluon(1) = 2
endif
if (k(lme1,2).eq.21) then
isgluon(2) = 1
elseif (abs(k(lme1,2)).le.5) then
isgluon(2) = 0
else
isgluon(2) = 2
endif
endif
ZA(LME1)=1.d0
ZA(LME2)=1.d0
THETAA(LME1)=P(LME1,5)/(SQRT(Z1*(1.-Z1))*P(LME1,4))
THETAA(LME2)=P(LME2,5)/(SQRT(Z2*(1.-Z2))*P(LME2,4))
ZD(LME1)=Z1
ZD(LME2)=Z2
QQBARD(LME1)=WHICH1
QQBARD(LME2)=WHICH2
MV(LME1,1)=X0
MV(LME1,2)=Y0
MV(LME1,3)=0.d0
MV(LME1,4)=0.d0
IF(P(LME1,5).GT.0.d0)THEN
LAMBDA=1.d0/(FTFAC*P(LME1,4)*0.2/Q1**2)
MV(LME1,5)=-LOG(1.d0-PYR(0))/LAMBDA
ELSE
MV(LME1,5)=LTIME
ENDIF
MV(LME2,1)=X0
MV(LME2,2)=Y0
MV(LME2,3)=0.d0
MV(LME2,4)=0.d0
IF(P(LME2,5).GT.0.d0)THEN
LAMBDA=1.d0/(FTFAC*P(LME2,4)*0.2/Q2**2)
MV(LME2,5)=-LOG(1.d0-PYR(0))/LAMBDA
ELSE
MV(LME2,5)=LTIME
ENDIF
C--develop parton shower
CALL MAKECASCADE
IF(DISCARD) THEN
NGOOD=NGOOD-1
WDISC=WDISC+EVWEIGHT
NDISC=NDISC+1
write(logfid,*)'discard event',J
GOTO 102
ENDIF
IF(.NOT.ALLHAD)THEN
DO 86 I=1,N
IF(K(I,1).EQ.3) K(I,1)=22
86 CONTINUE
ENDIF
IF(HADRO)THEN
CALL MAKESTRINGS(HADROTYPE)
IF(DISCARD) THEN
write(logfid,*)'discard event',J
WDISC=WDISC+EVWEIGHT
NDISC=NDISC+1
NGOOD=NGOOD-1
GOTO 102
ENDIF
CALL PYEXEC
IF(MSTU(30).NE.ERRCOUNT)THEN
write(logfid,*)'PYTHIA discards event',J,
& ' (error number',MSTU(30),')'
ERRCOUNT=MSTU(30)
WDISC=WDISC+EVWEIGHT
NDISC=NDISC+1
NGOOD=NGOOD-1
GOTO 102
ENDIF
ENDIF
DO 888 I=1,N
IF(K(I,2).EQ.94)THEN
NGOOD=NGOOD-1
NSTRANGE=NSTRANGE+1
NDISC=NDISC+1
GOTO 102
ENDIF
888 CONTINUE
IF(MSTU(30).NE.ERRCOUNT)THEN
ERRCOUNT=MSTU(30)
ELSE
CALL CONVERTTOHEPMC(HPMCFID,NGOOD,PID,b1,b2)
ENDIF
C--write message to log-file
102 IF(NSIM.GT.100)THEN
IF(MOD(J,NSIM/100).EQ.0)THEN
write(logfid,*) 'done with event number ',J,
& PARI(1), (sumofweights-wdisc)/j
ENDIF
else
write(logfid,*) 'done with event number ',J
ENDIF
call flush(logfid)
end
***********************************************************************
*** subroutine makestrings
***********************************************************************
SUBROUTINE MAKESTRINGS(WHICH)
IMPLICIT NONE
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
INTEGER WHICH
IF(WHICH.EQ.0)THEN
CALL MAKESTRINGS_VAC
ELSEIF(WHICH.EQ.1)THEN
CALL MAKESTRINGS_MINL
ELSE
WRITE(logfid,*)'error: unknown hadronisation type in MAKESTRINGS'
ENDIF
END
***********************************************************************
*** subroutine makestrings_vac
***********************************************************************
SUBROUTINE MAKESTRINGS_VAC
IMPLICIT NONE
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
C--Common block of Pythia
COMMON/PYJETS/N,NPAD,K(23000,5),P(23000,5),V(23000,5)
INTEGER N,NPAD,K
DOUBLE PRECISION P,V
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--colour index common block
COMMON/COLOUR/TRIP(23000),ANTI(23000),COLMAX
INTEGER TRIP,ANTI,COLMAX
C--discard event flag
COMMON/DISC/NDISC,NSTRANGE,NGOOD,errcount,wdisc,DISCARD
LOGICAL DISCARD
INTEGER NDISC,NSTRANGE,NGOOD,errcount
double precision wdisc
C--local variables
INTEGER NOLD,I,J,LQUARK,LMATCH,LLOOSE,NOLD1
DOUBLE PRECISION EADDEND,PYR,DIR
LOGICAL ISDIQUARK,compressevent,roomleft
DATA EADDEND/10.d0/
i = 0
if (compress) roomleft = compressevent(i)
NOLD1=N
C--remove all active lines that are leptons, gammas, hadrons etc.
DO 52 I=1,NOLD1
IF((K(I,1).EQ.4).AND.(TRIP(I).EQ.0).AND.(ANTI(I).EQ.0))THEN
C--copy line to end of event record
N=N+1
IF(N.GT.22990) THEN
write(logfid,*)'event too long for event record'
DISCARD=.TRUE.
RETURN
ENDIF
K(N,1)=11
K(N,2)=K(I,2)
K(N,3)=I
K(N,4)=0
K(N,5)=0
P(N,1)=P(I,1)
P(N,2)=P(I,2)
P(N,3)=P(I,3)
P(N,4)=P(I,4)
P(N,5)=P(I,5)
K(I,1)=17
K(I,4)=N
K(I,5)=N
TRIP(N)=TRIP(I)
ANTI(N)=ANTI(I)
ENDIF
52 CONTINUE
NOLD=N
C--first do strings with existing (anti)triplets
C--find string end (=quark or antiquark)
43 LQUARK=0
DO 40 I=1,NOLD
IF((K(I,1).EQ.11).OR.(K(I,1).EQ.12).OR.(K(I,1).EQ.13)
& .OR.(K(I,1).EQ.14)) K(I,1)=17
IF(((K(I,1).EQ.1).OR.(K(I,1).EQ.3).OR.(K(I,1).EQ.4))
& .AND.((K(I,2).LT.6).OR.ISDIQUARK(K(I,2))))THEN
LQUARK=I
GOTO 41
ENDIF
40 CONTINUE
GOTO 50
41 CONTINUE
C--copy string end to end of event record
N=N+1
IF(N.GT.22990) THEN
write(logfid,*)'event too long for event record'
DISCARD=.TRUE.
RETURN
ENDIF
K(N,1)=2
K(N,2)=K(LQUARK,2)
K(N,3)=LQUARK
K(N,4)=0
K(N,5)=0
P(N,1)=P(LQUARK,1)
P(N,2)=P(LQUARK,2)
P(N,3)=P(LQUARK,3)
P(N,4)=P(LQUARK,4)
P(N,5)=P(LQUARK,5)
K(LQUARK,1)=16
K(LQUARK,4)=N
K(LQUARK,5)=N
TRIP(N)=TRIP(LQUARK)
ANTI(N)=ANTI(LQUARK)
C--append matching colour partner
LMATCH=0
DO 44 J=1,10000000
DO 42 I=1,NOLD
IF(((K(I,1).EQ.1).OR.(K(I,1).EQ.3).OR.(K(I,1).EQ.4))
& .AND.(((TRIP(I).EQ.ANTI(N)).AND.(TRIP(I).NE.0))
& .OR.((ANTI(I).EQ.TRIP(N)).AND.(ANTI(I).NE.0))))THEN
N=N+1
IF(N.GT.22990) THEN
write(logfid,*)'event too long for event record'
DISCARD=.TRUE.
RETURN
ENDIF
K(N,2)=K(I,2)
K(N,3)=I
K(N,4)=0
K(N,5)=0
P(N,1)=P(I,1)
P(N,2)=P(I,2)
P(N,3)=P(I,3)
P(N,4)=P(I,4)
P(N,5)=P(I,5)
TRIP(N)=TRIP(I)
ANTI(N)=ANTI(I)
K(I,1)=16
K(I,4)=N
K(I,5)=N
IF(K(I,2).EQ.21)THEN
K(N,1)=2
GOTO 44
ELSE
K(N,1)=1
GOTO 43
ENDIF
ENDIF
42 CONTINUE
C--no matching colour partner found, add artificial end point
write(logfid,*)'Error in MAKESTRINGS_VAC: failed to reconstruct '//
&'colour singlet system, will discard event'
discard = .true.
return
44 CONTINUE
C--now take care of purely gluonic remainder system
C-----------------------------------------
C--find gluon where anti-triplet is not matched
50 LLOOSE=0
DO 45 I=1,NOLD
IF(((K(I,1).EQ.1).OR.(K(I,1).EQ.3).OR.(K(I,1).EQ.4)))THEN
DO 46 J=1,NOLD
IF(((K(I,1).EQ.1).OR.(K(I,1).EQ.3).OR.(K(I,1).EQ.4)))THEN
IF(ANTI(I).EQ.TRIP(J)) GOTO 45
ENDIF
46 CONTINUE
LLOOSE=I
GOTO 47
ENDIF
45 CONTINUE
GOTO 51
47 CONTINUE
C--generate artificial triplet end
write(logfid,*)'Error in MAKESTRINGS_VAC: failed to reconstruct '//
&'colour singlet system, will discard event'
discard = .true.
return
C--copy loose gluon to end of event record
N=N+1
IF(N.GT.22990) THEN
write(logfid,*)'event too long for event record'
DISCARD=.TRUE.
RETURN
ENDIF
K(N,1)=2
K(N,2)=K(LLOOSE,2)
K(N,3)=LLOOSE
K(N,4)=0
K(N,5)=0
P(N,1)=P(LLOOSE,1)
P(N,2)=P(LLOOSE,2)
P(N,3)=P(LLOOSE,3)
P(N,4)=P(LLOOSE,4)
P(N,5)=P(LLOOSE,5)
K(LLOOSE,1)=16
K(LLOOSE,4)=N
K(LLOOSE,5)=N
TRIP(N)=TRIP(LLOOSE)
ANTI(N)=ANTI(LLOOSE)
C--append matching colour partner
LMATCH=0
DO 48 J=1,10000000
DO 49 I=1,NOLD
IF(((K(I,1).EQ.1).OR.(K(I,1).EQ.3).OR.(K(I,1).EQ.4))
& .AND.(ANTI(I).EQ.TRIP(N)))THEN
N=N+1
IF(N.GT.22990) THEN
write(logfid,*)'event too long for event record'
DISCARD=.TRUE.
RETURN
ENDIF
K(N,2)=K(I,2)
K(N,3)=I
K(N,4)=0
K(N,5)=0
P(N,1)=P(I,1)
P(N,2)=P(I,2)
P(N,3)=P(I,3)
P(N,4)=P(I,4)
P(N,5)=P(I,5)
TRIP(N)=TRIP(I)
ANTI(N)=ANTI(I)
K(I,1)=16
K(I,4)=N
K(I,5)=N
K(N,1)=2
GOTO 48
ENDIF
49 CONTINUE
C--no matching colour partner found, add artificial end point
write(logfid,*)'Error in MAKESTRINGS_VAC: failed to reconstruct '//
&'colour singlet system, will discard event'
discard = .true.
return
48 CONTINUE
51 CONTINUE
CALL CLEANUP(NOLD1)
END
***********************************************************************
*** subroutine makestrings_minl
***********************************************************************
SUBROUTINE MAKESTRINGS_MINL
IMPLICIT NONE
C--Common block of Pythia
COMMON/PYJETS/N,NPAD,K(23000,5),P(23000,5),V(23000,5)
INTEGER N,NPAD,K
DOUBLE PRECISION P,V
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--local variables
INTEGER NOLD,I,J,LMAX,LMIN,LEND
DOUBLE PRECISION EMAX,MINV,MMIN,Z,GENERATEZ,MCUT,EADDEND,PYR,DIR
DATA MCUT/1.d8/
DATA EADDEND/10.d0/
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
C--discard event flag
COMMON/DISC/NDISC,NSTRANGE,NGOOD,errcount,wdisc,DISCARD
LOGICAL DISCARD
INTEGER NDISC,NSTRANGE,NGOOD,errcount
double precision wdisc
logical compressevent,roomleft
i = 0
if (compress) roomleft = compressevent(i)
NOLD=N
C--find most energetic unfragmented parton in event
43 EMAX=0
LMAX=0
DO 40 I=1,NOLD
IF((K(I,1).EQ.11).OR.(K(I,1).EQ.12).OR.(K(I,1).EQ.13)
& .OR.(K(I,1).EQ.14)) K(I,1)=17
IF(((K(I,1).EQ.1).OR.(K(I,1).EQ.3).OR.(K(I,1).EQ.4))
& .AND.(P(I,4).GT.EMAX))THEN
EMAX=P(I,4)
LMAX=I
ENDIF
40 CONTINUE
C--if there is non, we are done
IF(LMAX.EQ.0) GOTO 50
C--check if highest energy parton is (anti)quark or gluon
IF(K(LMAX,2).EQ.21)THEN
C--split gluon in qqbar pair and store one temporarily in line 1
C--make new line in event record for string end
N=N+2
IF(N.GT.22990) THEN
write(logfid,*)'event too long for event record'
DISCARD=.TRUE.
RETURN
ENDIF
IF((N-2).GT.NOLD)THEN
DO 47 J=NOLD,N-3
K(N+NOLD-J,1)=K(N+NOLD-J-2,1)
K(N+NOLD-J,2)=K(N+NOLD-J-2,2)
IF(K(N+NOLD-J-2,3).GT.NOLD) THEN
K(N+NOLD-J,3)=K(N+NOLD-J-2,3)+2
ELSE
K(N+NOLD-J,3)=K(N+NOLD-J-2,3)
ENDIF
K(N+NOLD-J,4)=0
K(N+NOLD-J,5)=0
P(N+NOLD-J,1)=P(N+NOLD-J-2,1)
P(N+NOLD-J,2)=P(N+NOLD-J-2,2)
P(N+NOLD-J,3)=P(N+NOLD-J-2,3)
P(N+NOLD-J,4)=P(N+NOLD-J-2,4)
P(N+NOLD-J,5)=P(N+NOLD-J-2,5)
K(K(N+NOLD-J-2,3),4)=K(K(N+NOLD-J-2,3),4)+2
K(K(N+NOLD-J-2,3),5)=K(K(N+NOLD-J-2,3),5)+2
47 CONTINUE
ENDIF
NOLD=NOLD+2
K(LMAX,1)=18
Z=GENERATEZ(0.d0,0.d0,1.d-3,'QG')
IF(Z.GT.0.5)THEN
K(NOLD-1,2)=1
K(NOLD,2)=-1
ELSE
Z=1.-Z
K(NOLD-1,2)=-1
K(NOLD,2)=1
ENDIF
K(NOLD-1,1)=1
K(NOLD-1,3)=LMAX
K(NOLD-1,4)=0
K(NOLD-1,5)=0
P(NOLD-1,1)=(1.-Z)*P(LMAX,1)
P(NOLD-1,2)=(1.-Z)*P(LMAX,2)
P(NOLD-1,3)=(1.-Z)*P(LMAX,3)
P(NOLD-1,4)=(1.-Z)*P(LMAX,4)
P(NOLD-1,5)=P(LMAX,5)
K(NOLD,1)=1
K(NOLD,3)=LMAX
K(NOLD,4)=0
K(NOLD,5)=0
P(NOLD,1)=Z*P(LMAX,1)
P(NOLD,2)=Z*P(LMAX,2)
P(NOLD,3)=Z*P(LMAX,3)
P(NOLD,4)=Z*P(LMAX,4)
P(NOLD,5)=P(LMAX,5)
K(LMAX,1)=18
K(LMAX,4)=NOLD-1
K(LMAX,5)=NOLD
LMAX=NOLD
ENDIF
N=N+1
IF(N.GT.22990) THEN
write(logfid,*)'event too long for event record'
DISCARD=.TRUE.
RETURN
ENDIF
K(N,1)=2
K(N,2)=K(LMAX,2)
K(N,3)=LMAX
K(N,4)=0
K(N,5)=0
P(N,1)=P(LMAX,1)
P(N,2)=P(LMAX,2)
P(N,3)=P(LMAX,3)
P(N,4)=P(LMAX,4)
P(N,5)=P(LMAX,5)
K(LMAX,1)=16
K(LMAX,4)=N
K(LMAX,5)=N
LEND=LMAX
C--find closest partner
42 MMIN=1.d10
LMIN=0
DO 41 I=1,NOLD
IF(((K(I,1).EQ.1).OR.(K(I,1).EQ.3).OR.(K(I,1).EQ.4))
& .AND.((K(I,2).EQ.21).OR.((K(I,2)*K(LEND,2).LT.0.d0).AND.
& (K(I,3).NE.K(LEND,3))))
& .AND.(P(I,1)*P(LEND,1).GT.0.d0))THEN
MINV=P(I,4)*P(LMAX,4)-P(I,1)*P(LMAX,1)-P(I,2)*P(LMAX,2)
& -P(I,3)*P(LMAX,3)
IF((MINV.LT.MMIN).AND.(MINV.GT.0.d0).AND.(MINV.LT.MCUT))THEN
MMIN=MINV
LMIN=I
ENDIF
ENDIF
41 CONTINUE
C--if no closest partner can be found, generate artificial end point for string
IF(LMIN.EQ.0)THEN
N=N+1
IF(N.GT.22990) THEN
write(logfid,*)'event too long for event record'
DISCARD=.TRUE.
RETURN
ENDIF
K(N,1)=1
K(N,2)=-K(LEND,2)
K(N,3)=0
K(N,4)=0
K(N,5)=0
P(N,1)=0.d0
P(N,2)=0.d0
IF(PYR(0).LT.0.5)THEN
DIR=1.d0
ELSE
DIR=-1.d0
ENDIF
P(N,3)=DIR*EADDEND
P(N,4)=EADDEND
P(N,5)=0.d0
GOTO 43
ELSE
C--else build closest partner in string
N=N+1
IF(N.GT.22990) THEN
write(logfid,*)'event too long for event record'
DISCARD=.TRUE.
RETURN
ENDIF
K(N,2)=K(LMIN,2)
K(N,3)=LMIN
K(N,4)=0
K(N,5)=0
P(N,1)=P(LMIN,1)
P(N,2)=P(LMIN,2)
P(N,3)=P(LMIN,3)
P(N,4)=P(LMIN,4)
P(N,5)=P(LMIN,5)
K(LMIN,1)=16
K(LMIN,4)=N
K(LMIN,5)=N
IF(K(LMIN,2).EQ.21)THEN
K(N,1)=2
LMAX=LMIN
GOTO 42
ELSE
K(N,1)=1
GOTO 43
ENDIF
ENDIF
50 CONTINUE
CALL CLEANUP(NOLD)
END
***********************************************************************
*** subroutine cleanup
***********************************************************************
SUBROUTINE CLEANUP(NFIRST)
IMPLICIT NONE
C--Common block of Pythia
COMMON/PYJETS/N,NPAD,K(23000,5),P(23000,5),V(23000,5)
INTEGER N,NPAD,K
DOUBLE PRECISION P,V
C--local variables
INTEGER NFIRST,NLAST,I,J
NLAST=N
DO 21 I=1,NLAST-NFIRST
DO 22 J=1,5
K(I,J)=K(NFIRST+I,J)
P(I,J)=P(NFIRST+I,J)
V(I,J)=V(NFIRST+I,J)
22 CONTINUE
K(I,3)=0
21 CONTINUE
N=NLAST-NFIRST
END
***********************************************************************
*** subroutine makecascade
***********************************************************************
SUBROUTINE MAKECASCADE
IMPLICIT NONE
C--Common block of Pythia
COMMON/PYJETS/N,NPAD,K(23000,5),P(23000,5),V(23000,5)
INTEGER N,NPAD,K
DOUBLE PRECISION P,V
C--time common block
COMMON/TIME/MV(23000,5)
DOUBLE PRECISION MV
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--discard event flag
COMMON/DISC/NDISC,NSTRANGE,NGOOD,errcount,wdisc,DISCARD
LOGICAL DISCARD
INTEGER NDISC,NSTRANGE,NGOOD,errcount
double precision wdisc
C--local variables
INTEGER NOLD,I
LOGICAL CONT
10 NOLD=N
CONT=.FALSE.
DO 11 I=2,NOLD
if (i.gt.n) goto 10
C--check if parton may evolve, i.e. do splitting or scattering
IF((K(I,1).EQ.1).OR.(K(I,1).EQ.2))THEN
CONT=.TRUE.
CALL MAKEBRANCH(I)
IF(DISCARD) GOTO 12
ENDIF
11 CONTINUE
IF(CONT) GOTO 10
12 END
***********************************************************************
*** subroutine makebranch
***********************************************************************
SUBROUTINE MAKEBRANCH(L)
IMPLICIT NONE
C--Common block of Pythia
COMMON/PYJETS/N,NPAD,K(23000,5),P(23000,5),V(23000,5)
INTEGER N,NPAD,K
DOUBLE PRECISION P,V
C--time common block
COMMON/TIME/MV(23000,5)
DOUBLE PRECISION MV
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--discard event flag
COMMON/DISC/NDISC,NSTRANGE,NGOOD,errcount,wdisc,DISCARD
LOGICAL DISCARD
INTEGER NDISC,NSTRANGE,NGOOD,errcount
double precision wdisc
C--variables for angular ordering
COMMON/ANGOR/ZA(23000),ZD(23000),THETAA(23000),QQBARD(23000)
DOUBLE PRECISION ZA,ZD,THETAA
LOGICAL QQBARD
C--number of scattering events
COMMON/CHECK/NSCAT,NSCATEFF,NSPLIT,nspliti,nsplitf,nistry,
&nisfail,nfsfail,nfstry,nttot,ntrej
DOUBLE PRECISION NSCAT,NSCATEFF,NSPLIT,nspliti,nsplitf,nistry,
&nisfail,nfsfail,nfstry,nttot,ntrej
C--variables for coherent scattering
COMMON/COHERENT/NSTART,NEND,ALLQS(10000,6),SCATCENTRES(10000,10),
&QSUMVEC(4),QSUM2
INTEGER NSTART,NEND
DOUBLE PRECISION ALLQS,SCATCENTRES,QSUMVEC,QSUM2
C--event weight
COMMON/WEIGHT/EVWEIGHT,sumofweights
double precision EVWEIGHT,sumofweights
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
C--local variables
INTEGER L,LINE,NOLD,TYPI,LINEOLD,LKINE,nendold
DOUBLE PRECISION THETA,PHI,PYP,FORMTIME,STARTTIME,TLEFT,
&TSUM,DELTAT,NEWMASS,GETMASS,Q,GETMS,ZDEC,X,MS,DTCORR
LOGICAL OVERQ0,QQBARDEC
CHARACTER TYP
LOGICAL RADIATION,RETRYSPLIT,MEDIND,roomleft,compressevent
LINE=L
NSTART=0
NEND=0
STARTTIME=MV(LINE,4)
TSUM=0.d0
QSUM2=0.d0
QSUMVEC(1)=0.d0
QSUMVEC(2)=0.d0
QSUMVEC(3)=0.d0
QSUMVEC(4)=0.d0
RETRYSPLIT=.FALSE.
MEDIND=.FALSE.
X=0.d0
Q=0.d0
TYPI=0
IF ((N.GT.20000).and.compress) roomleft = compressevent(line)
20 IF(DISCARD) RETURN
MS=GETMS(MV(LINE,1),MV(LINE,2),MV(LINE,3),MV(LINE,4))
IF(((K(LINE,1).EQ.1).AND.(P(LINE,5).GT.0.d0))
& .OR.((K(LINE,1).EQ.2).AND.(P(LINE,5).GT.MS)))THEN
IF(MEDIND)THEN
FORMTIME=starttime
ELSE
FORMTIME=MIN(MV(LINE,5),LTIME)
ENDIF
RADIATION=.TRUE.
ELSE
FORMTIME=LTIME
RADIATION=.FALSE.
ENDIF
TLEFT=FORMTIME-STARTTIME
IF(K(LINE,2).EQ.21)THEN
TYP='G'
ELSE
TYP='Q'
ENDIF
MEDIND=.FALSE.
IF(TLEFT.LE.1.d-10)THEN
C--no scattering
IF(RADIATION)THEN
C--if there is radiation associated with the parton then form it now
C--rotate such that momentum points in z-direction
NOLD=N
THETA=PYP(LINE,13)
PHI=PYP(LINE,15)
CALL PYROBO(LINE,LINE,0d0,-PHI,0d0,0d0,0d0)
CALL PYROBO(LINE,LINE,-THETA,0d0,0d0,0d0,0d0)
CALL MAKESPLITTING(LINE)
C--rotate back
CALL PYROBO(LINE,LINE,THETA,0d0,0d0,0d0,0d0)
CALL PYROBO(LINE,LINE,0d0,PHI,0d0,0d0,0d0)
IF(DISCARD) RETURN
CALL PYROBO(N-1,N,THETA,0d0,0d0,0d0,0d0)
CALL PYROBO(N-1,N,0d0,PHI,0d0,0d0,0d0)
C--set the production vertices: x_mother + (tprod - tprod_mother) * beta_mother
MV(N-1,1)=MV(LINE,1)
& +(MV(N-1,4)-MV(LINE,4))*P(LINE,1)/max(pyp(line,8),P(LINE,4))
MV(N-1,2)=MV(LINE,2)
& +(MV(N-1,4)-MV(LINE,4))*P(LINE,2)/max(pyp(line,8),P(LINE,4))
MV(N-1,3)=MV(LINE,3)
& +(MV(N-1,4)-MV(LINE,4))*P(LINE,3)/max(pyp(line,8),P(LINE,4))
MV(N, 1)=MV(LINE,1)
& +(MV(N, 4)-MV(LINE,4))*P(LINE,1)/max(pyp(line,8),P(LINE,4))
MV(N, 2)=MV(LINE,2)
& +(MV(N, 4)-MV(LINE,4))*P(LINE,2)/max(pyp(line,8),P(LINE,4))
MV(N, 3)=MV(LINE,3)
& +(MV(N, 4)-MV(LINE,4))*P(LINE,3)/max(pyp(line,8),P(LINE,4))
LINE=N
NSTART=0
NEND=0
STARTTIME=MV(N,4)
QSUMVEC(1)=0.d0
QSUMVEC(2)=0.d0
QSUMVEC(3)=0.d0
QSUMVEC(4)=0.d0
QSUM2=0.d0
TSUM=0.d0
GOTO 21
ELSE
NSTART=0
NEND=0
STARTTIME=FORMTIME
QSUMVEC(1)=0.d0
QSUMVEC(2)=0.d0
QSUMVEC(3)=0.d0
QSUMVEC(4)=0.d0
QSUM2=0.d0
TSUM=0.d0
GOTO 21
ENDIF
ELSE
C--do scattering
C--find delta t for the scattering
DELTAT=TLEFT
OVERQ0=.FALSE.
CALL DOINSTATESCAT(LINE,X,TYPI,Q,STARTTIME+TSUM,DELTAT,
& OVERQ0,.FALSE.)
TSUM=TSUM+DELTAT
TLEFT=TLEFT-DELTAT
C--do initial state splitting if there is one
NOLD=N
LINEOLD=LINE
ZDEC=ZD(LINE)
QQBARDEC=QQBARD(LINE)
25 IF(X.LT.1.d0) THEN
CALL MAKEINSPLIT(LINE,X,QSUM2,Q,TYPI,STARTTIME+TSUM,DELTAT)
IF(DISCARD) RETURN
IF(X.LT.1.d0)THEN
LINE=N
LKINE=N
IF(K(LINE,2).EQ.21)THEN
NEWMASS=GETMASS(0.d0,SCALEFACM*SQRT(-QSUM2),-1.d0,P(LINE,4),
& 21,SQRT(-QSUM2),.FALSE.,ZDEC,QQBARDEC)
IF(ZDEC.GT.0.d0)THEN
THETAA(LINE)=NEWMASS/(SQRT(ZDEC*(1.-ZDEC))*P(LINE,4))
ELSE
THETAA(LINE)=0.d0
ENDIF
ZD(LINE)=ZDEC
QQBARD(LINE)=QQBARDEC
ELSE
NEWMASS=GETMASS(0.d0,SCALEFACM*SQRT(-QSUM2),-1.d0,P(LINE,4),
& k(line,2),SQRT(-QSUM2),.FALSE.,ZDEC,QQBARDEC)
IF(ZDEC.GT.0.d0)THEN
THETAA(LINE)=NEWMASS/(SQRT(ZDEC*(1.-ZDEC))*P(LINE,4))
ELSE
THETAA(LINE)=0.d0
ENDIF
ZD(LINE)=ZDEC
QQBARD(LINE)=QQBARDEC
ENDIF
ZDEC=ZD(LINE)
QQBARDEC=QQBARD(LINE)
ELSE
LKINE=LINE
NEND=NSTART
QSUM2=ALLQS(NEND,1)
QSUMVEC(1)=ALLQS(NEND,2)
QSUMVEC(2)=ALLQS(NEND,3)
QSUMVEC(3)=ALLQS(NEND,4)
QSUMVEC(4)=ALLQS(NEND,5)
IF(-ALLQS(NEND,1).GT.Q0**2/SCALEFACM**2)THEN
OVERQ0=.TRUE.
ELSE
OVERQ0=.FALSE.
ENDIF
tleft = starttime+tsum+tleft-allqs(1,6)
tsum = allqs(1,6)-starttime
ENDIF
ENDIF
IF(X.EQ.1.d0)THEN
NEWMASS=0.d0
IF(NEND.GT.0)THEN
CALL DOFISTATESCAT(LINE,STARTTIME+TSUM,TLEFT,DELTAT,
& NEWMASS,OVERQ0,ZDEC,QQBARDEC)
IF(NEWMASS.GT.(P(LINE,5)+1.d-10))THEN
MEDIND=.TRUE.
ELSE
MEDIND=.FALSE.
ZDEC=ZD(LINE)
QQBARDEC=QQBARD(LINE)
ENDIF
TSUM=TSUM+DELTAT
TLEFT=TLEFT-DELTAT
LKINE=LINE
ENDIF
ENDIF
C--do kinematics
RETRYSPLIT=.FALSE.
IF(NEND.GT.0) THEN
nendold=nend
CALL DOKINEMATICS(LKINE,lineold,NSTART,NEND,NEWMASS,RETRYSPLIT,
& STARTTIME+TSUM,X,ZDEC,QQBARDEC)
IF(RETRYSPLIT) THEN
tleft = starttime+tsum+tleft-allqs(1,6)
tsum = allqs(1,6)-starttime
if (x.lt.1.d0) then
NEND=NSTART
QSUM2=ALLQS(NEND,1)
QSUMVEC(1)=ALLQS(NEND,2)
QSUMVEC(2)=ALLQS(NEND,3)
QSUMVEC(3)=ALLQS(NEND,4)
QSUMVEC(4)=ALLQS(NEND,5)
TYPI=K(L,2)
IF(-ALLQS(NEND,1).GT.Q0**2/SCALEFACM**2)THEN
OVERQ0=.TRUE.
ELSE
OVERQ0=.FALSE.
ENDIF
N=NOLD
LINE=LINEOLD
X=1.d0
K(LINE,1)=1
NSPLIT=NSPLIT-EVWEIGHT
nspliti=nspliti-evweight
GOTO 25
else
LINE=N
STARTTIME=STARTTIME+TSUM
TSUM=0.d0
endif
ELSE
LINE=N
STARTTIME=STARTTIME+TSUM
TSUM=0.d0
ENDIF
ELSE
STARTTIME=STARTTIME+TSUM
TSUM=0.d0
ENDIF
IF(P(LINE,5).GT.0.d0) RADIATION=.TRUE.
ENDIF
21 IF(((K(LINE,1).EQ.1).AND.(P(LINE,5).GT.0.d0))
& .OR.((K(LINE,1).EQ.2).AND.(P(LINE,5).GT.MS))
& .OR.(STARTTIME.LT.LTIME))THEN
GOTO 20
ENDIF
IF((K(LINE,1).EQ.1).AND.(P(LINE,5).EQ.0.d0)) K(LINE,1)=4
IF((K(LINE,1).EQ.2).AND.(P(LINE,5).LE.MS)) K(LINE,1)=4
END
***********************************************************************
*** subroutine makesplitting
***********************************************************************
SUBROUTINE MAKESPLITTING(L)
IMPLICIT NONE
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
C--Common block of Pythia
COMMON/PYJETS/N,NPAD,K(23000,5),P(23000,5),V(23000,5)
INTEGER N,NPAD,K
DOUBLE PRECISION P,V
C--time common block
COMMON/TIME/MV(23000,5)
DOUBLE PRECISION MV
C--factor in front of formation times
COMMON/FTIMEFAC/FTFAC
DOUBLE PRECISION FTFAC
C--colour index common block
COMMON/COLOUR/TRIP(23000),ANTI(23000),COLMAX
INTEGER TRIP,ANTI,COLMAX
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--discard event flag
COMMON/DISC/NDISC,NSTRANGE,NGOOD,errcount,wdisc,DISCARD
LOGICAL DISCARD
INTEGER NDISC,NSTRANGE,NGOOD,errcount
double precision wdisc
C--variables for angular ordering
COMMON/ANGOR/ZA(23000),ZD(23000),THETAA(23000),QQBARD(23000)
DOUBLE PRECISION ZA,ZD,THETAA
LOGICAL QQBARD
C--number of scattering events
COMMON/CHECK/NSCAT,NSCATEFF,NSPLIT,nspliti,nsplitf,nistry,
&nisfail,nfsfail,nfstry,nttot,ntrej
DOUBLE PRECISION NSCAT,NSCATEFF,NSPLIT,nspliti,nsplitf,nistry,
&nisfail,nfsfail,nfstry,nttot,ntrej
C--event weight
COMMON/WEIGHT/EVWEIGHT,sumofweights
double precision EVWEIGHT,sumofweights
C--local variables
INTEGER L,DIR,quarkflav
DOUBLE PRECISION PHIQ,PYR,PI,GENERATEZ,BMAX1,CMAX1,PTS,MB,MC,
&GETMASS,PZ,EPS,QH,Z,R,LAMBDA,WEIGHT,ZDECB,ZDECC,XDEC(3),THETA,
&GETTEMP
LOGICAL QUARK,QQBAR,QQBARDECB,QQBARDECC
integer bin
DATA PI/3.141592653589793d0/
IF((N+2).GT.22990) THEN
write(logfid,*)'event too long for event record'
DISCARD=.TRUE.
RETURN
ENDIF
XDEC(1)=MV(L,1)+(MV(L,5)-MV(L,4))*P(L,1)/P(L,4)
XDEC(2)=MV(L,2)+(MV(L,5)-MV(L,4))*P(L,2)/P(L,4)
XDEC(3)=MV(L,3)+(MV(L,5)-MV(L,4))*P(L,3)/P(L,4)
IF(GETTEMP(XDEC(1),XDEC(2),XDEC(3),MV(L,5)).GT.0.d0)THEN
THETA=-1.d0
ELSE
THETA=THETAA(L)
ENDIF
C--on-shell partons cannot split
IF((P(L,5).EQ.0d0).OR.(K(L,1).EQ.11).OR.(K(L,1).EQ.12)
& .OR.(K(L,1).EQ.13).OR.(K(L,1).EQ.14).OR.(K(L,1).EQ.3)) GOTO 31
C--quark or gluon?
IF(K(L,2).EQ.21)THEN
QUARK=.FALSE.
ELSE
QUARK=.TRUE.
QQBAR=.FALSE.
quarkflav=k(L,2)
ENDIF
C--if gluon decide on kind of splitting
QQBAR=QQBARD(L)
if (qqbar) quarkflav=int(pyr(0)*3.)+1
C--if g->gg splitting decide on colour order
IF(QUARK.OR.QQBAR)THEN
DIR=0
ELSE
IF(PYR(0).LT.0.5)THEN
DIR=1
ELSE
DIR=-1
ENDIF
ENDIF
Z=ZD(L)
IF(Z.EQ.0.d0)THEN
write(logfid,*)'makesplitting: z=0',L
goto 36
ENDIF
GOTO 35
C--generate z value
36 IF(ANGORD.AND.(ZA(L).NE.1.d0))THEN
C--additional z constraint due to angular ordering
QH=4.*P(L,5)**2*(1.-ZA(L))/(ZA(L)*P(K(L,3),5)**2)
IF(QH.GT.1)THEN
write(logfid,*)L,': reject event: angular ordering
& conflict in medium'
CALL PYLIST(3)
DISCARD=.TRUE.
GOTO 31
ENDIF
EPS=0.5-0.5*SQRT(1.-QH)
ELSE
EPS=0d0
ENDIF
IF(QUARK)THEN
Z=GENERATEZ(P(L,5)**2,P(L,4),EPS,'QQ')
ELSE
IF(QQBAR)THEN
Z=GENERATEZ(P(L,5)**2,P(L,4),EPS,'QG')
ELSE
Z=GENERATEZ(P(L,5)**2,P(L,4),EPS,'GG')
ENDIF
ENDIF
35 CONTINUE
C--maximum virtualities for daughters
BMAX1=MIN(P(L,5),Z*P(L,4))
CMAX1=MIN(P(L,5),(1.-Z)*P(L,4))
C--generate mass of quark or gluon (particle b) from Sudakov FF
30 IF(QUARK.OR.QQBAR)THEN
MB=GETMASS(0.d0,BMAX1,THETA,Z*P(L,4),quarkflav,
& BMAX1,.FALSE.,ZDECB,QQBARDECB)
ELSE
MB=GETMASS(0.d0,BMAX1,THETA,Z*P(L,4),21,
& BMAX1,.FALSE.,ZDECB,QQBARDECB)
ENDIF
C--generate mass gluon (particle c) from Sudakov FF
IF(QUARK.OR.(.NOT.QQBAR))THEN
MC=GETMASS(0.d0,CMAX1,THETA,(1.-Z)*P(L,4),21,
& CMAX1,.FALSE.,ZDECC,QQBARDECC)
ELSE
MC=GETMASS(0.d0,CMAX1,THETA,(1.-Z)*P(L,4),-quarkflav,
& CMAX1,.FALSE.,ZDECC,QQBARDECC)
ENDIF
C--quark (parton b) momentum
182 PZ=(2.*Z*P(L,4)**2-P(L,5)**2-MB**2+MC**2)/(2.*P(L,3))
PTS=Z**2*(P(L,4)**2)-PZ**2-MB**2
C--if kinematics doesn't work out, generate new virtualities
C for daughters
C--massive phase space weight
IF((MB.EQ.0.d0).AND.(MC.EQ.0.d0).AND.(PTS.LT.0.d0)) GOTO 36
WEIGHT=1.d0
IF((PYR(0).GT.WEIGHT).OR.(PTS.LT.0.d0)
& .OR.((MB+MC).GT.P(L,5)))THEN
IF(MB.GT.MC)THEN
IF(QUARK.OR.QQBAR)THEN
MB=GETMASS(0.d0,MB,THETA,Z*P(L,4),quarkflav,
& BMAX1,.FALSE.,ZDECB,QQBARDECB)
ELSE
MB=GETMASS(0.d0,MB,THETA,Z*P(L,4),21,
& BMAX1,.FALSE.,ZDECB,QQBARDECB)
ENDIF
ELSE
IF(QUARK.OR.(.NOT.QQBAR))THEN
MC=GETMASS(0.d0,MC,THETA,(1.-Z)*P(L,4),21,
& CMAX1,.FALSE.,ZDECC,QQBARDECC)
ELSE
MC=GETMASS(0.d0,MC,THETA,(1.-Z)*P(L,4),-quarkflav,
& CMAX1,.FALSE.,ZDECC,QQBARDECC)
ENDIF
ENDIF
GOTO 182
ENDIF
N=N+2
C--take care of first daughter (radiated gluon or antiquark)
K(N-1,1)=K(L,1)
IF(QQBAR)THEN
K(N-1,2)=-quarkflav
TRIP(N-1)=0
ANTI(N-1)=ANTI(L)
ELSE
K(N-1,2)=21
IF((K(L,2).GT.0).AND.(DIR.GE.0))THEN
TRIP(N-1)=TRIP(L)
ANTI(N-1)=COLMAX+1
ELSE
TRIP(N-1)=COLMAX+1
ANTI(N-1)=ANTI(L)
ENDIF
COLMAX=COLMAX+1
ENDIF
K(N-1,3)=L
K(N-1,4)=0
K(N-1,5)=0
P(N-1,4)=(1-Z)*P(L,4)
P(N-1,5)=MC
ZA(N-1)=1.-Z
IF(ZDECC.GT.0.d0)THEN
THETAA(N-1)=P(N-1,5)/(SQRT(ZDECC*(1.-ZDECC))*P(N-1,4))
ELSE
THETAA(N-1)=0.d0
ENDIF
ZD(N-1)=ZDECC
QQBARD(N-1)=QQBARDECC
C--take care of second daughter (final quark or gluon or quark from
C gluon splitting)
K(N,1)=K(L,1)
IF(QUARK)THEN
K(N,2)=K(L,2)
IF(K(N,2).GT.0)THEN
TRIP(N)=ANTI(N-1)
ANTI(N)=0
ELSE
TRIP(N)=0
ANTI(N)=TRIP(N-1)
ENDIF
ELSEIF(QQBAR)THEN
K(N,2)=quarkflav
TRIP(N)=TRIP(L)
ANTI(N)=0
ELSE
K(N,2)=21
IF(DIR.EQ.1)THEN
TRIP(N)=ANTI(N-1)
ANTI(N)=ANTI(L)
ELSE
TRIP(N)=TRIP(L)
ANTI(N)=TRIP(N-1)
ENDIF
ENDIF
K(N,3)=L
K(N,4)=0
K(N,5)=0
P(N,3)=PZ
P(N,4)=Z*P(L,4)
P(N,5)=MB
ZA(N)=Z
IF(ZDECB.GT.0.d0)THEN
THETAA(N)=P(N,5)/(SQRT(ZDECB*(1.-ZDECB))*P(N,4))
ELSE
THETAA(N)=0.d0
ENDIF
ZD(N)=ZDECB
QQBARD(N)=QQBARDECB
C--azimuthal angle
PHIQ=2*PI*PYR(0)
P(N,1)=SQRT(PTS)*COS(PHIQ)
P(N,2)=SQRT(PTS)*SIN(PHIQ)
C--gluon momentum
P(N-1,1)=P(L,1)-P(N,1)
P(N-1,2)=P(L,2)-P(N,2)
P(N-1,3)=P(L,3)-P(N,3)
MV(N-1,4)=MV(L,5)
IF(P(N-1,5).GT.0.d0)THEN
LAMBDA=1.d0/(FTFAC*P(N-1,4)*0.2/P(N-1,5)**2)
MV(N-1,5)=MV(L,5)-LOG(1.d0-PYR(0))/LAMBDA
ELSE
MV(N-1,5)=0.d0
ENDIF
MV(N,4)=MV(L,5)
IF(P(N,5).GT.0.d0)THEN
LAMBDA=1.d0/(FTFAC*P(N,4)*0.2/P(N,5)**2)
MV(N,5)=MV(L,5)-LOG(1.d0-PYR(0))/LAMBDA
ELSE
MV(N,5)=0.d0
ENDIF
C--take care of initial quark (or gluon)
IF(K(L,1).EQ.2)THEN
K(L,1)=13
ELSE
K(L,1)=11
ENDIF
K(L,4)=N-1
K(L,5)=N
NSPLIT=NSPLIT+EVWEIGHT
nsplitf=nsplitf+evweight
31 CONTINUE
END
***********************************************************************
*** subroutine makeinsplit
***********************************************************************
SUBROUTINE MAKEINSPLIT(L,X,TSUM,VIRT,TYPI,TIME,TAURAD)
IMPLICIT NONE
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
C--Common block of Pythia
COMMON/PYJETS/N,NPAD,K(23000,5),P(23000,5),V(23000,5)
INTEGER N,NPAD,K
DOUBLE PRECISION P,V
C--time common block
COMMON/TIME/MV(23000,5)
DOUBLE PRECISION MV
C--factor in front of formation times
COMMON/FTIMEFAC/FTFAC
DOUBLE PRECISION FTFAC
C--colour index common block
COMMON/COLOUR/TRIP(23000),ANTI(23000),COLMAX
INTEGER TRIP,ANTI,COLMAX
C--variables for angular ordering
COMMON/ANGOR/ZA(23000),ZD(23000),THETAA(23000),QQBARD(23000)
DOUBLE PRECISION ZA,ZD,THETAA
LOGICAL QQBARD
C--discard event flag
COMMON/DISC/NDISC,NSTRANGE,NGOOD,errcount,wdisc,DISCARD
LOGICAL DISCARD
INTEGER NDISC,NSTRANGE,NGOOD,errcount
double precision wdisc
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--number of scattering events
COMMON/CHECK/NSCAT,NSCATEFF,NSPLIT,nspliti,nsplitf,nistry,
&nisfail,nfsfail,nfstry,nttot,ntrej
DOUBLE PRECISION NSCAT,NSCATEFF,NSPLIT,nspliti,nsplitf,nistry,
&nisfail,nfsfail,nfstry,nttot,ntrej
C--event weight
COMMON/WEIGHT/EVWEIGHT,sumofweights
double precision EVWEIGHT,sumofweights
C--local variables
INTEGER L,TYPI,NOLD,DIR,typc
DOUBLE PRECISION X,VIRT,MB2,MC2,GETMASS,PZ,KT2,THETA,PHI,PI,
&PHIQ,PYP,PYR,R,TIME,TSUM,TAURAD,LAMBDA,ZDEC
LOGICAL QQBARDEC
CHARACTER*2 TYP2
integer bin
DATA PI/3.141592653589793d0/
IF((N+2).GT.22990) THEN
write(logfid,*)'event too long for event record'
DISCARD=.TRUE.
RETURN
ENDIF
IF(K(L,2).EQ.21)THEN
IF(TYPI.EQ.21)THEN
TYP2='GG'
TYPC=21
ELSE
TYP2='QG'
TYPC=typi
ENDIF
ELSE
IF(TYPI.EQ.21)THEN
TYP2='GQ'
TYPC=k(L,2)
ELSE
TYP2='QQ'
TYPC=21
ENDIF
ENDIF
C--if g->gg decide on colour configuration
IF(TYP2.EQ.'GG')THEN
IF(PYR(0).LT.0.5)THEN
DIR=1
ELSE
DIR=-1
ENDIF
ELSE
DIR=0
ENDIF
MB2=VIRT**2
MB2=P(L,5)**2-MB2
MC2=GETMASS(0.d0,SCALEFACM*SQRT(-TSUM),-1.d0,
& (1.-X)*P(L,4),TYPC,(1.-X)*P(L,4),
& .FALSE.,ZDEC,QQBARDEC)**2
C--rotate such that momentum points in z-direction
NOLD=N
THETA=PYP(L,13)
PHI=PYP(L,15)
CALL PYROBO(L,L,0d0,-PHI,0d0,0d0,0d0)
CALL PYROBO(L,L,-THETA,0d0,0d0,0d0,0d0)
PZ=(2*X*P(L,4)**2-P(L,5)**2-MB2+MC2)/(2*P(L,3))
KT2=X**2*(P(L,4)**2)-PZ**2-MB2
IF(KT2.LT.0.d0)THEN
MC2=0.d0
PZ=(2*X*P(L,4)**2-P(L,5)**2-MB2+MC2)/(2*P(L,3))
KT2=X**2*(P(L,4)**2)-PZ**2-MB2
IF(KT2.LT.0.d0)THEN
CALL PYROBO(L,L,THETA,0d0,0d0,0d0,0d0)
CALL PYROBO(L,L,0d0,PHI,0d0,0d0,0d0)
X=1.d0
RETURN
ENDIF
ENDIF
N=N+2
C--take care of first daughter (radiated gluon or antiquark)
K(N-1,1)=K(L,1)
IF(TYP2.EQ.'QG')THEN
K(N-1,2)=-TYPI
IF(K(N-1,2).GT.0)THEN
TRIP(N-1)=TRIP(L)
ANTI(N-1)=0
ELSE
TRIP(N-1)=0
ANTI(N-1)=ANTI(L)
ENDIF
ELSEIF(TYP2.EQ.'GQ')THEN
K(N-1,2)=K(L,2)
IF(K(N-1,2).GT.0)THEN
TRIP(N-1)=COLMAX+1
ANTI(N-1)=0
ELSE
TRIP(N-1)=0
ANTI(N-1)=COLMAX+1
ENDIF
COLMAX=COLMAX+1
ELSE
K(N-1,2)=21
IF((K(L,2).GT.0).AND.(DIR.GE.0))THEN
TRIP(N-1)=TRIP(L)
ANTI(N-1)=COLMAX+1
ELSE
TRIP(N-1)=COLMAX+1
ANTI(N-1)=ANTI(L)
ENDIF
COLMAX=COLMAX+1
ENDIF
K(N-1,3)=L
K(N-1,4)=0
K(N-1,5)=0
P(N-1,4)=(1.-X)*P(L,4)
P(N-1,5)=SQRT(MC2)
C--take care of second daughter (final quark or gluon or quark from
C gluon splitting)
K(N,1)=K(L,1)
IF(TYP2.EQ.'QG')THEN
K(N,2)=TYPI
IF(K(N,2).GT.0)THEN
TRIP(N)=TRIP(L)
ANTI(N)=0
ELSE
TRIP(N)=0
ANTI(N)=ANTI(L)
ENDIF
ELSEIF(TYPI.NE.21)THEN
K(N,2)=K(L,2)
IF(K(N,2).GT.0)THEN
TRIP(N)=ANTI(N-1)
ANTI(N)=0
ELSE
TRIP(N)=0
ANTI(N)=TRIP(N-1)
ENDIF
ELSE
K(N,2)=21
IF(K(N-1,2).EQ.21)THEN
IF(DIR.EQ.1)THEN
TRIP(N)=ANTI(N-1)
ANTI(N)=ANTI(L)
ELSE
TRIP(N)=TRIP(L)
ANTI(N)=TRIP(N-1)
ENDIF
ELSEIF(K(N-1,2).GT.0)THEN
TRIP(N)=TRIP(L)
ANTI(N)=TRIP(N-1)
ELSE
TRIP(N)=ANTI(N-1)
ANTI(N)=ANTI(L)
ENDIF
ENDIF
K(N,3)=L
K(N,4)=0
K(N,5)=0
P(N,3)=PZ
P(N,4)=X*P(L,4)
IF(MB2.LT.0.d0)THEN
P(N,5)=-SQRT(-MB2)
ELSE
P(N,5)=SQRT(MB2)
ENDIF
C--azimuthal angle
PHIQ=2*PI*PYR(0)
P(N,1)=SQRT(KT2)*COS(PHIQ)
P(N,2)=SQRT(KT2)*SIN(PHIQ)
C--gluon momentum
P(N-1,1)=P(L,1)-P(N,1)
P(N-1,2)=P(L,2)-P(N,2)
P(N-1,3)=P(L,3)-P(N,3)
MV(L,5)=TIME-TAURAD
MV(N-1,4)=MV(L,5)
IF(P(N-1,5).GT.0.d0)THEN
LAMBDA=1.d0/(FTFAC*P(N-1,4)*0.2/P(N-1,5)**2)
MV(N-1,5)=MV(L,5)-LOG(1.d0-PYR(0))/LAMBDA
ELSE
MV(N-1,5)=0.d0
ENDIF
MV(N,4)=MV(L,5)
IF(P(N,5).GT.0.d0)THEN
MV(N,5)=TIME
ELSE
MV(N,5)=0.d0
ENDIF
ZA(N-1)=1.d0
THETAA(N-1)=-1.d0
ZD(N-1)=ZDEC
QQBARD(N-1)=QQBARDEC
ZA(N)=1.d0
THETAA(N)=-1.d0
ZD(N)=0.d0
QQBARD(N)=.FALSE.
C--take care of initial quark (or gluon)
IF(K(L,1).EQ.2)THEN
K(L,1)=13
ELSE
K(L,1)=11
ENDIF
K(L,4)=N-1
K(L,5)=N
NSPLIT=NSPLIT+EVWEIGHT
nspliti=nspliti+evweight
CALL PYROBO(L,L,THETA,0d0,0d0,0d0,0d0)
CALL PYROBO(N-1,N,THETA,0d0,0d0,0d0,0d0)
CALL PYROBO(L,L,0d0,PHI,0d0,0d0,0d0)
CALL PYROBO(N-1,N,0d0,PHI,0d0,0d0,0d0)
C--set the production vertices: x_mother + (tprod - tprod_mother) * beta_mother
MV(N-1,1)=MV(L,1)+(MV(N-1,4)-MV(L,4))*P(L,1)/max(pyp(l,8),P(L,4))
MV(N-1,2)=MV(L,2)+(MV(N-1,4)-MV(L,4))*P(L,2)/max(pyp(l,8),P(L,4))
MV(N-1,3)=MV(L,3)+(MV(N-1,4)-MV(L,4))*P(L,3)/max(pyp(l,8),P(L,4))
MV(N, 1)=MV(L,1)+(MV(N, 4)-MV(L,4))*P(L,1)/max(pyp(l,8),P(L,4))
MV(N, 2)=MV(L,2)+(MV(N, 4)-MV(L,4))*P(L,2)/max(pyp(l,8),P(L,4))
MV(N, 3)=MV(L,3)+(MV(N, 4)-MV(L,4))*P(L,3)/max(pyp(l,8),P(L,4))
END
***********************************************************************
*** subroutine doinstatescat
***********************************************************************
SUBROUTINE DOINSTATESCAT(L,X,TYPI,Q,TSTART,DELTAT,OVERQ0,
& RETRYSPLIT)
IMPLICIT NONE
C--Common block of Pythia
COMMON/PYJETS/N,NPAD,K(23000,5),P(23000,5),V(23000,5)
INTEGER N,NPAD,K
DOUBLE PRECISION P,V
C--time common block
COMMON/TIME/MV(23000,5)
DOUBLE PRECISION MV
C--factor in front of formation times
COMMON/FTIMEFAC/FTFAC
DOUBLE PRECISION FTFAC
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--discard event flag
COMMON/DISC/NDISC,NSTRANGE,NGOOD,errcount,wdisc,DISCARD
LOGICAL DISCARD
INTEGER NDISC,NSTRANGE,NGOOD,errcount
double precision wdisc
C--variables for coherent scattering
COMMON/COHERENT/NSTART,NEND,ALLQS(10000,6),SCATCENTRES(10000,10),
&QSUMVEC(4),QSUM2
INTEGER NSTART,NEND
DOUBLE PRECISION ALLQS,SCATCENTRES,QSUMVEC,QSUM2
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
C--local variables
INTEGER L,TYPI,COUNTER,COUNTMAX,COUNT2
DOUBLE PRECISION X,DELTAT,DELTAL,PYR,R,PNORAD,GETPNORAD1,GETNOSCAT,
&WEIGHT,LOW,FMAX,GETPDF,SIGMATOT,GETSSCAT,PFCHANGE,PI,TNOW,TLEFT,
&XMAX,PQQ,PQG,PGQ,PGG,ALPHAS,TSTART,TSUM,Q,QOLD,Q2OLD,GETNEWMASS,
&GENERATEZ,TMAX,TMAXNEW,DT,XSC,YSC,ZSC,TSC,MS1,MD1,GETMS,GETMD,
&GETTEMP,GETNEFF,LAMBDA,RTAU,PHI,TAUEST,QSUMVECOLD(4),ZDUM,WEIGHT,
&pyp
LOGICAL FCHANGE,NORAD,OVERQ0,NOSCAT,GETDELTAT,RETRYSPLIT,
&QQBARDUM
CHARACTER TYP
CHARACTER*2 TYP2
DATA PI/3.141592653589793d0/
DATA COUNTMAX/10000/
COUNTER=0
XSC=MV(L,1)+(TSTART-MV(L,4))*P(L,1)/P(L,4)
YSC=MV(L,2)+(TSTART-MV(L,4))*P(L,2)/P(L,4)
ZSC=MV(L,3)+(TSTART-MV(L,4))*P(L,3)/P(L,4)
TSC=TSTART
MD1=GETMD(XSC,YSC,ZSC,TSC)
MS1=GETMS(XSC,YSC,ZSC,TSC)
IF(MD1.LE.1.D-4.OR.MS1.LE.1.D-4)THEN
write(logfid,*)'problem!',GETTEMP(XSC,YSC,ZSC,TSC),
&GETNEFF(XSC,YSC,ZSC,TSC)
ENDIF
C--check for scattering
NOSCAT=.NOT.GETDELTAT(L,TSTART,DELTAT,DT)
IF(NOSCAT.AND.(.NOT.RETRYSPLIT)) GOTO 116
C--decide whether there will be radiation
PNORAD=GETPNORAD1(L,xsc,ysc,zsc,tsc)
IF((PYR(0).LT.PNORAD).OR.(P(L,4).LT.1.001*Q0))THEN
NORAD=.TRUE.
ELSE
NORAD=.FALSE.
ENDIF
C--decide whether q or g is to be scattered
IF(K(L,2).EQ.21)THEN
TYP='G'
TYP2='GC'
SIGMATOT=GETSSCAT(P(L,4),p(l,1),p(l,2),p(l,3),P(L,5),
& Q0,'G','C',xsc,ysc,zsc,tsc,0)
IF((SIGMATOT.EQ.0.d0).OR.(PNORAD.EQ.1.d0))THEN
PFCHANGE=0.d0
ELSE
PFCHANGE=GETSSCAT(P(L,4),p(l,1),p(l,2),p(l,3),P(L,5),
& Q0,'G','Q',xsc,ysc,zsc,tsc,0)
& /SIGMATOT
ENDIF
SIGMATOT=GETSSCAT(P(L,4),p(l,1),p(l,2),p(l,3),P(L,5),
& 0.d0,'G','C',xsc,ysc,zsc,tsc,0)
ELSE
TYP='Q'
TYP2='QQ'
SIGMATOT=GETSSCAT(P(L,4),p(l,1),p(l,2),p(l,3),P(L,5),
& Q0,'Q','C',xsc,ysc,zsc,tsc,0)
IF((SIGMATOT.EQ.0.d0).OR.(PNORAD.EQ.1.d0))THEN
PFCHANGE=0.d0
ELSE
PFCHANGE=GETSSCAT(P(L,4),p(l,1),p(l,2),p(l,3),P(L,5),
& Q0,'Q','G',xsc,ysc,zsc,tsc,0)
& /SIGMATOT
ENDIF
SIGMATOT=GETSSCAT(P(L,4),p(l,1),p(l,2),p(l,3),P(L,5),
& 0.d0,'Q','C',xsc,ysc,zsc,tsc,0)
ENDIF
IF((PFCHANGE.LT.-1.d-4).OR.(PFCHANGE.GT.1.d0+1.d-4)) THEN
write(logfid,*)'error: flavour change probability=',
& PFCHANGE,'for ',TYP
ENDIF
IF(PYR(0).LT.PFCHANGE)THEN
FCHANGE=.TRUE.
ELSE
FCHANGE=.FALSE.
ENDIF
IF (NORAD) FCHANGE=.FALSE.
C--set TYPI
IF(TYP.EQ.'G')THEN
IF(FCHANGE)THEN
typi = int(pyr(0)*3)+1
r=pyr(0)
if (r.gt.0.5) typi=-typi
ELSE
TYPI=K(L,2)
ENDIF
ELSE
IF(FCHANGE)THEN
TYPI=21
ELSE
TYPI=K(L,2)
ENDIF
ENDIF
LOW=Q0**2/SCALEFACM**2
TMAX=4.*(P(L,4)**2-P(L,5)**2)
XMAX=1.-Q0**2/(SCALEFACM**2*4.*TMAX)
IF(SIGMATOT.EQ.0.d0) GOTO 116
RTAU=PYR(0)
C--generate a trial emission
C--pick a x value from splitting function
112 COUNTER=COUNTER+1
IF(TYP.EQ.'G')THEN
IF(FCHANGE)THEN
X=GENERATEZ(0.d0,0.d0,1.-XMAX,'QG')
ELSE
X=GENERATEZ(0.d0,0.d0,1.-XMAX,'GG')
ENDIF
ELSE
IF(FCHANGE)THEN
X=1.-GENERATEZ(0.d0,0.d0,1.-XMAX,'QQ')
ELSE
X=GENERATEZ(0.d0,0.d0,1.-XMAX,'QQ')
ENDIF
ENDIF
IF(NORAD) X=1.d0
C--initialisation
TMAXNEW=(X*P(L,4))**2
PHI=0.d0
TLEFT=DELTAT
TNOW=TSTART
QSUMVEC(1)=0.d0
QSUMVEC(2)=0.d0
QSUMVEC(3)=0.d0
QSUMVEC(4)=0.d0
QSUM2=-1.d-10
OVERQ0=.FALSE.
Q=P(L,5)
QOLD=P(L,5)
TAUEST=DELTAT
C--generate first momentum transfer
DELTAL=DT
NSTART=1
NEND=1
TNOW=TNOW+DELTAL
TSUM=DELTAL
TLEFT=TLEFT-DELTAL
ALLQS(NEND,6)=TNOW
Q2OLD=QSUM2
C--get new momentum transfer
COUNT2=0
118 CALL GETQVEC(L,NEND,TNOW-MV(L,4),X)
IF(-QSUM2.GT.P(L,4)**2)THEN
QSUMVEC(1)=0.d0
QSUMVEC(2)=0.d0
QSUMVEC(3)=0.d0
QSUMVEC(4)=0.d0
QSUM2=Q2OLD
IF(COUNT2.LT.100)THEN
COUNT2=COUNT2+1
GOTO 118
ELSE
ALLQS(NEND,1)=0.d0
ALLQS(NEND,2)=0.d0
ALLQS(NEND,3)=0.d0
ALLQS(NEND,4)=0.d0
ALLQS(NEND,5)=0.d0
ENDIF
ENDIF
C--update OVERQ0
IF(-ALLQS(NEND,1).GT.LOW) OVERQ0=.TRUE.
C--get new virtuality
IF(OVERQ0.AND.(.NOT.NORAD))THEN
Q=GETNEWMASS(L,SCALEFACM**2*QSUM2,SCALEFACM**2*Q2OLD,0.d0,
& .TRUE.,X,ZDUM,QQBARDUM)
ELSE
Q=0.d0
ENDIF
C--estimate formation time
111 IF((Q.EQ.0.d0).OR.(Q.EQ.P(L,5)))THEN
TAUEST=DELTAT
ELSE
TAUEST=FTFAC*(1.-PHI)*0.2*X*P(L,4)/Q**2
ENDIF
LAMBDA=1.d0/TAUEST
TAUEST=-LOG(1.d0-RTAU)/LAMBDA
C--find number, position and momentum transfers of further scatterings
NOSCAT=.NOT.GETDELTAT(L,TNOW,MIN(TLEFT,TAUEST),DELTAL)
IF((.NOT.NOSCAT).AND.(.NOT.NORAD))THEN
C--add a momentum transfer
NEND=NEND+1
IF(NEND.GE.100)THEN
nend=nend-1
goto 114
ENDIF
TNOW=TNOW+DELTAL
TSUM=TSUM+DELTAL
TLEFT=TLEFT-DELTAL
C--update phase
IF((Q.NE.0.d0).AND.(Q.NE.P(L,5)))THEN
PHI=PHI+5.*DELTAL*Q**2/(1.*X*P(L,4))
ENDIF
C--get new momentum transfer
ALLQS(NEND,6)=TNOW
Q2OLD=QSUM2
QSUMVECOLD(1)=QSUMVEC(1)
QSUMVECOLD(2)=QSUMVEC(2)
QSUMVECOLD(3)=QSUMVEC(3)
QSUMVECOLD(4)=QSUMVEC(4)
COUNT2=0
119 CALL GETQVEC(L,NEND,TNOW-MV(L,4),X)
IF(-QSUM2.GT.P(L,4)**2)THEN
QSUMVEC(1)=QSUMVECOLD(1)
QSUMVEC(2)=QSUMVECOLD(2)
QSUMVEC(3)=QSUMVECOLD(3)
QSUMVEC(4)=QSUMVECOLD(4)
QSUM2=Q2OLD
IF(COUNT2.LT.100)THEN
COUNT2=COUNT2+1
GOTO 119
ELSE
ALLQS(NEND,1)=0.d0
ALLQS(NEND,2)=0.d0
ALLQS(NEND,3)=0.d0
ALLQS(NEND,4)=0.d0
ALLQS(NEND,5)=0.d0
ENDIF
ENDIF
C--update OVERQ0
IF((-QSUM2.GT.LOW)
& .OR.(-ALLQS(NEND,1).GT.LOW)) OVERQ0=.TRUE.
C--get new virtuality
QOLD=Q
IF(OVERQ0.AND.(.NOT.NORAD))THEN
Q=GETNEWMASS(L,SCALEFACM**2*QSUM2,SCALEFACM**2*Q2OLD,0.d0,
& .TRUE.,X,ZDUM,QQBARDUM)
ELSE
Q=0.d0
ENDIF
GOTO 111
ENDIF
C--do reweighting
114 TMAXNEW=X**2*P(L,4)**2
IF(NORAD)THEN
WEIGHT=1.d0
Q=0.d0
X=1.d0
ELSEIF((-QSUM2.LT.LOW).OR.(Q.EQ.0.d0))THEN
WEIGHT=0.d0
ELSEIF(-QSUM2.GT.P(L,4)**2)THEN
WEIGHT=0.d0
ELSE
IF(TYP.EQ.'G')THEN
FMAX=2.*LOG(-SCALEFACM**2*QSUM2/Q0**2)
& *ALPHAS(Q0**2/4.,LPS)/(2.*PI)
IF(QSUM2.EQ.0.d0)THEN
WEIGHT=0.d0
NORAD=.TRUE.
ELSE
IF(FCHANGE)THEN
WEIGHT=2.*GETPDF(X,SCALEFACM*SQRT(-QSUM2),'QG')/(PQG(X)*FMAX)
IF((WEIGHT.GT.1.d0+1.d-4).OR.(WEIGHT.LT.-1.d-4))THEN
write(logfid,*)'x,sqrt(qsum^2),getpdf,fmax:',X,
& SQRT(-QSUM2),GETPDF(X,SCALEFACM*SQRT(-QSUM2),'QG'),'qg',
& FMAX
ENDIF
ELSE
WEIGHT=GETPDF(X,SCALEFACM*SQRT(-QSUM2),'GG')/(PGG(X)*FMAX)
IF((WEIGHT.GT.1.d0+1.d-4).OR.(WEIGHT.LT.-1.d-4))THEN
write(logfid,*)'x,sqrt(qsum^2),getpdf,fmax:',X,
& SQRT(-QSUM2),GETPDF(X,SCALEFACM*SQRT(-QSUM2),'GG'),'gg',
& FMAX
ENDIF
ENDIF
ENDIF
ELSE
FMAX=LOG(-SCALEFACM**2*QSUM2/Q0**2)
& *ALPHAS(Q0**2/4.,LPS)/(2.*PI)
IF(QSUM2.EQ.0.d0)THEN
WEIGHT=0.d0
NORAD=.TRUE.
ELSE
IF(FCHANGE)THEN
WEIGHT=GETPDF(X,SCALEFACM*SQRT(-QSUM2),'GQ')/(PGQ(X)*FMAX)
IF((WEIGHT.GT.1.d0+1.d-4).OR.(WEIGHT.LT.-1.d-4))THEN
write(logfid,*)'x,sqrt(qsum^2),getpdf:,fmax',X,
& SQRT(-QSUM2),GETPDF(X,SCALEFACM*SQRT(-QSUM2),'GQ'),'gq',
& FMAX
ENDIF
ELSE
WEIGHT=GETPDF(X,SCALEFACM*SQRT(-QSUM2),'QQ')/(PQQ(X)*FMAX)
IF((WEIGHT.GT.1.d0+1.d-4).OR.(WEIGHT.LT.-1.d-4))THEN
write(logfid,*)'x,sqrt(qsum^2),getpdf,fmax:',X,
& SQRT(-QSUM2),GETPDF(X,SCALEFACM*SQRT(-QSUM2),'QQ'),'qq',
& FMAX
ENDIF
ENDIF
ENDIF
ENDIF
ENDIF
IF((WEIGHT.GT.1.d0+1.d-4).OR.(WEIGHT.LT.-1.d-4))
& write(logfid,*)'error: weight=',WEIGHT
115 IF(PYR(0).GT.WEIGHT)THEN
IF(COUNTER.LT.COUNTMAX)THEN
GOTO 112
ELSE
Q=0.d0
X=1.d0
NEND=NSTART
QSUM2=ALLQS(NEND,1)
QSUMVEC(1)=ALLQS(NEND,2)
QSUMVEC(2)=ALLQS(NEND,3)
QSUMVEC(3)=ALLQS(NEND,4)
QSUMVEC(4)=ALLQS(NEND,5)
TYPI=K(L,2)
IF(-ALLQS(NEND,1).GT.LOW)THEN
OVERQ0=.TRUE.
ELSE
OVERQ0=.FALSE.
ENDIF
DELTAT=ALLQS(NEND,6)-TSTART
TNOW=ALLQS(1,6)
RETURN
ENDIF
ENDIF
C--found meaningful configuration, now do final checks
C--check if phase is unity and weight with 1/Nscat
IF(((TLEFT.LT.TAUEST).OR.(PYR(0).GT.1.d0/(NEND*1.d0)))
& .AND.(.NOT.NORAD))THEN
Q=0.d0
X=1.d0
NEND=NSTART
QSUM2=ALLQS(NEND,1)
QSUMVEC(1)=ALLQS(NEND,2)
QSUMVEC(2)=ALLQS(NEND,3)
QSUMVEC(3)=ALLQS(NEND,4)
QSUMVEC(4)=ALLQS(NEND,5)
TYPI=K(L,2)
IF(-ALLQS(NEND,1).GT.LOW)THEN
OVERQ0=.TRUE.
ELSE
OVERQ0=.FALSE.
ENDIF
DELTAT=ALLQS(NEND,6)-TSTART
TNOW=ALLQS(1,6)
ELSE
IF(.NOT.NORAD)THEN
TLEFT=TLEFT-TAUEST
TNOW=TNOW+TAUEST
TSUM=TSUM+TAUEST
ENDIF
DELTAT=TSUM
ENDIF
RETURN
C--exit in case of failure
116 Q=0.d0
X=1.d0
NSTART=0
NEND=0
QSUMVEC(1)=0.d0
QSUMVEC(2)=0.d0
QSUMVEC(3)=0.d0
QSUMVEC(4)=0.d0
QSUM2=0.d0
OVERQ0=.FALSE.
TYPI=K(L,2)
RETURN
END
***********************************************************************
*** subroutine dofistatescat
***********************************************************************
SUBROUTINE DOFISTATESCAT(L,TNOW,DTLEFT,DELTAT,NEWMASS,
& OVERQ0,Z,QQBAR)
IMPLICIT NONE
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
C--Common block of Pythia
COMMON/PYJETS/N,NPAD,K(23000,5),P(23000,5),V(23000,5)
INTEGER N,NPAD,K
DOUBLE PRECISION P,V
C--time common block
COMMON/TIME/MV(23000,5)
DOUBLE PRECISION MV
C--factor in front of formation times
COMMON/FTIMEFAC/FTFAC
DOUBLE PRECISION FTFAC
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--discard event flag
COMMON/DISC/NDISC,NSTRANGE,NGOOD,errcount,wdisc,DISCARD
LOGICAL DISCARD
INTEGER NDISC,NSTRANGE,NGOOD,errcount
double precision wdisc
C--variables for coherent scattering
COMMON/COHERENT/NSTART,NEND,ALLQS(10000,6),SCATCENTRES(10000,10),
&QSUMVEC(4),QSUM2
INTEGER NSTART,NEND
DOUBLE PRECISION ALLQS,SCATCENTRES,QSUMVEC,QSUM2
C--local variables
INTEGER L,COUNTER,COUNTMAX,COUNT2
DOUBLE PRECISION TNOW,DELTAT,NEWMASS,TLEFT,DELTAL,Q2OLD,
&GETNEWMASS,PYR,TSUM,QSUMVECOLD(4),RTAU,LAMBDA,DTLEFT,PHI,
&TAUEST,LOW,Z,pyp
LOGICAL OVERQ0,NOSCAT,GETDELTAT,QQBAR
CHARACTER TYP
DATA COUNTMAX/100/
DELTAL=0.d0
IF(-QSUM2.GT.P(L,4)**2)
& write(logfid,*) 'DOFISTATESCAT has a problem:',-QSUM2,P(L,4)**2
IF(K(L,2).EQ.21)THEN
TYP='G'
ELSE
TYP='Q'
ENDIF
LOW=Q0**2/SCALEFACM**2
TSUM=0.d0
PHI=0.d0
DELTAT=0.d0
C--check for radiation with first (given) momentum transfer
Q2OLD=0.d0
IF(OVERQ0.OR.(-QSUM2.GT.LOW))THEN
NEWMASS=GETNEWMASS(L,SCALEFACM**2*QSUM2,SCALEFACM**2*Q2OLD,
& NEWMASS,.FALSE.,1.d0,Z,QQBAR)
OVERQ0=.TRUE.
ELSE
NEWMASS=P(L,5)
ENDIF
RTAU=PYR(0)
TLEFT=DTLEFT
222 IF((NEWMASS.EQ.0.d0).OR.(NEWMASS.EQ.P(L,5)))THEN
TAUEST=TLEFT
ELSE
TAUEST=FTFAC*(1.-PHI)*0.2*P(L,4)/NEWMASS**2
ENDIF
LAMBDA=1.d0/TAUEST
TAUEST=-LOG(1.d0-RTAU)/LAMBDA
NOSCAT=.NOT.GETDELTAT(L,TNOW+TSUM,MIN(TAUEST,TLEFT),DELTAL)
IF(.NOT.NOSCAT)THEN
C--do scattering
NEND=NEND+1
IF(NEND.gt.countmax)THEN
nend=nend-1
goto 218
ENDIF
IF(NSTART.EQ.0) NSTART=1
TSUM=TSUM+DELTAL
TLEFT=TLEFT-DELTAL
IF((NEWMASS.NE.0.d0).AND.(NEWMASS.NE.P(L,5)))THEN
PHI=PHI+5.*DELTAL*NEWMASS**2/(1.*P(L,4))
ENDIF
ALLQS(NEND,6)=TNOW+TSUM
QSUMVECOLD(1)=QSUMVEC(1)
QSUMVECOLD(2)=QSUMVEC(2)
QSUMVECOLD(3)=QSUMVEC(3)
QSUMVECOLD(4)=QSUMVEC(4)
Q2OLD=QSUM2
C--get new momentum transfer
COUNT2=0
219 CALL GETQVEC(L,NEND,TNOW+TSUM-MV(L,4),1.d0)
IF(-QSUM2.GT.P(L,4)**2)THEN
QSUMVEC(1)=QSUMVECOLD(1)
QSUMVEC(2)=QSUMVECOLD(2)
QSUMVEC(3)=QSUMVECOLD(3)
QSUMVEC(4)=QSUMVECOLD(4)
QSUM2=Q2OLD
IF(COUNT2.LT.100)THEN
COUNT2=COUNT2+1
GOTO 219
ELSE
ALLQS(NEND,1)=0.d0
ALLQS(NEND,2)=0.d0
ALLQS(NEND,3)=0.d0
ALLQS(NEND,4)=0.d0
ALLQS(NEND,5)=0.d0
ENDIF
ENDIF
C--figure out new virtuality
IF(OVERQ0.OR.(-QSUM2.GT.LOW))THEN
NEWMASS=GETNEWMASS(L,SCALEFACM**2*QSUM2,SCALEFACM**2*Q2OLD,
& NEWMASS,.FALSE.,1.d0,Z,QQBAR)
OVERQ0=.TRUE.
ENDIF
GOTO 222
ENDIF
C--no more scattering
218 if ((newmass**2.gt.low).and.(newmass.ne.p(l,5))) then
if ((TLEFT.LT.TAUEST).OR.(PYR(0).GT.1.d0/(NEND*1.d0))) then
if (nend.eq.countmax) then
deltat=tsum
else if (TLEFT.LT.TAUEST) then
DELTAT=TSUM+tleft
else
DELTAT=TSUM+tauest
endif
NEWMASS=P(L,5)
ELSE
DELTAT=TSUM+TAUEST
ENDIF
else
DELTAT=0.d0
NSTART=1
NEND=1
QSUM2=ALLQS(NEND,1)
QSUMVEC(1)=ALLQS(NEND,2)
QSUMVEC(2)=ALLQS(NEND,3)
QSUMVEC(3)=ALLQS(NEND,4)
QSUMVEC(4)=ALLQS(NEND,5)
IF(-ALLQS(NEND,1).GT.LOW)THEN
OVERQ0=.TRUE.
ELSE
OVERQ0=.FALSE.
ENDIF
NEWMASS=P(L,5)
endif
return
END
***********************************************************************
*** function getnewmass
***********************************************************************
DOUBLE PRECISION FUNCTION GETNEWMASS(L,Q2,QOLD2,MASS,IN,X,
& ZDEC,QQBARDEC)
IMPLICIT NONE
C--Common block of Pythia
COMMON/PYJETS/N,NPAD,K(23000,5),P(23000,5),V(23000,5)
INTEGER N,NPAD,K
DOUBLE PRECISION P,V
C--time common block
COMMON/TIME/MV(23000,5)
DOUBLE PRECISION MV
C--variables for angular ordering
COMMON/ANGOR/ZA(23000),ZD(23000),THETAA(23000),QQBARD(23000)
DOUBLE PRECISION ZA,ZD,THETAA
LOGICAL QQBARD
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--local variables
INTEGER L,flav
DOUBLE PRECISION Q2,QOLD2,R,PYR,PNOSPLIT1,PNOSPLIT2,Z,QA,
&GETSUDAKOV,GETMASS,PKEEP,X,MASS,ZDEC,QTMP,ZOLD
LOGICAL IN,QQBARDEC,QQBAROLD
CHARACTER*2 TYP
IF(x*P(L,4).LT.Q0)THEN
GETNEWMASS=0.d0
ZDEC=0.d0
QQBARDEC=.FALSE.
RETURN
ENDIF
IF (-Q2.LT.Q0**2)THEN
GETNEWMASS=0.d0
RETURN
ENDIF
flav = k(l,2)
IF(K(L,2).EQ.21)THEN
TYP='GC'
ELSE
TYP='QQ'
ENDIF
IF(SQRT(-QOLD2).LE.Q0)THEN
IF(IN)THEN
GETNEWMASS=GETMASS(0.d0,SQRT(-Q2),-1.d0,
& X*P(L,4),flav,X*P(L,4),IN,ZDEC,QQBARDEC)
ELSE
GETNEWMASS=GETMASS(0.d0,SQRT(-Q2),-1.d0,P(L,4),flav,
& SQRT(-Q2),IN,ZDEC,QQBARDEC)
ENDIF
GETNEWMASS=MIN(GETNEWMASS,X*P(L,4))
RETURN
ENDIF
Z=1.d0
QA=1.d0
IF(MAX(P(L,5),MASS).GT.0.d0)THEN
IF(-Q2.GT.-QOLD2)THEN
ZOLD=ZDEC
QQBAROLD=QQBARDEC
- QTMP=GETMASS(SQRT(-QOLD2),SQRT(-Q2),-1.d0,X*P(L,4),flav,
+ QTMP=GETMASS(0.d0,SQRT(-Q2),-1.d0,X*P(L,4),flav,
& SQRT(-Q2),IN,ZDEC,QQBARDEC)
- IF(QTMP.EQ.0.d0)THEN
+ IF(QTMP.LT.SQRT(-QOLD2))THEN
GETNEWMASS=MASS
ZDEC=ZOLD
QQBARDEC=QQBAROLD
ELSE
GETNEWMASS=QTMP
ENDIF
ELSE
PNOSPLIT1=GETSUDAKOV(SQRT(-QOLD2),QA,Q0,Z,X*P(L,4),
& TYP,MV(L,4),IN)
PNOSPLIT2=GETSUDAKOV(SQRT(-Q2),QA,Q0,Z,X*P(L,4),
& TYP,MV(L,4),IN)
PKEEP=(1.-PNOSPLIT2)/(1.-PNOSPLIT1)
IF(PYR(0).LT.PKEEP)THEN
IF(P(L,5).LT.SQRT(-Q2))THEN
GETNEWMASS=MASS
ELSE
- 55 GETNEWMASS=GETMASS(0.d0,SQRT(-Q2),-1.d0,X*P(L,4),flav,
+ 55 GETNEWMASS=GETMASS(Q0,SQRT(-Q2),-1.d0,X*P(L,4),flav,
& SQRT(-Q2),IN,ZDEC,QQBARDEC)
IF((GETNEWMASS.EQ.0.d0).AND.(X*P(L,4).GT.Q0)) GOTO 55
ENDIF
ELSE
GETNEWMASS=0.d0
ZDEC=0.d0
QQBARDEC=.FALSE.
ENDIF
ENDIF
ELSE
IF(-Q2.GT.-QOLD2)THEN
- GETNEWMASS=GETMASS(MAX(SQRT(-QOLD2),Q0),SQRT(-Q2),-1.d0,
+ GETNEWMASS=GETMASS(0.d0,SQRT(-Q2),-1.d0,
& X*P(L,4),flav,X*P(L,4),IN,ZDEC,QQBARDEC)
+ if(getnewmass.lt.SQRT(-QOLD2))then
+ GETNEWMASS=0.d0
+ ZDEC=0.d0
+ QQBARDEC=.FALSE.
+ endif
ELSE
GETNEWMASS=0.d0
ZDEC=0.d0
QQBARDEC=.FALSE.
ENDIF
ENDIF
GETNEWMASS=MIN(GETNEWMASS,x*P(L,4))
END
***********************************************************************
*** function getpnorad1
***********************************************************************
DOUBLE PRECISION FUNCTION GETPNORAD1(LINE,x,y,z,t)
IMPLICIT NONE
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
C--Common block of Pythia
COMMON/PYJETS/N,NPAD,K(23000,5),P(23000,5),V(23000,5)
INTEGER N,NPAD,K
DOUBLE PRECISION P,V
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--local variables
INTEGER LINE
DOUBLE PRECISION UP,LOW,CCOL,SIGMATOT,GETSSCAT,GETXSECINT,
&SCATPRIMFUNC,MS1,MD1,shat,pcms2,avmom(5),x,y,z,t,getmd
md1 = getmd(x,y,z,t)
call avscatcen(x,y,z,t,
&avmom(1),avmom(2),avmom(3),avmom(4),avmom(5))
ms1 = avmom(5)
shat = avmom(5)**2 + p(line,5)**2 + 2.*(avmom(4)*p(line,4)
& -avmom(1)*p(line,1)-avmom(2)*p(line,2)-avmom(3)*p(line,3))
pcms2 = (shat+p(line,5)**2-ms1**2)**2/(4.*shat)-p(line,5)**2
up = 4.*pcms2
LOW=Q0**2/SCALEFACM**2
IF((UP.LE.LOW).OR.(P(LINE,4).LT.Q0/SCALEFACM))THEN
GETPNORAD1=1.d0
RETURN
ENDIF
IF(K(LINE,2).EQ.21)THEN
CCOL=3./2.
C--probability for no initial state radiation
SIGMATOT=GETSSCAT(P(LINE,4),p(line,1),p(line,2),p(line,3),
& P(LINE,5),0.d0,'G','C',x,y,z,t,0)
IF(SIGMATOT.EQ.0.d0)THEN
GETPNORAD1=-1.d0
RETURN
ENDIF
GETPNORAD1=(CCOL*(SCATPRIMFUNC(LOW,MD1)-
&SCATPRIMFUNC(0.d0,MD1))
& + GETXSECINT(UP,MD1,'GB'))/SIGMATOT
ELSE
CCOL=2./3.
C--probability for no initial state radiation
SIGMATOT=GETSSCAT(P(LINE,4),p(line,1),p(line,2),p(line,3),
& P(LINE,5),0.d0,'Q','C',x,y,z,t,0)
IF(SIGMATOT.EQ.0.d0)THEN
GETPNORAD1=1.d0
RETURN
ENDIF
GETPNORAD1=(CCOL*(SCATPRIMFUNC(LOW,MD1)-
&SCATPRIMFUNC(0.d0,MD1))
& + GETXSECINT(UP,MD1,'QB'))/SIGMATOT
ENDIF
IF((GETPNORAD1.LT.-1.d-4).OR.(GETPNORAD1.GT.1.d0+1.d-4))THEN
write(logfid,*)'error: P_norad=',GETPNORAD1,
& P(LINE,4),P(LINE,5),LOW,UP,K(LINE,2),MD1
ENDIF
END
***********************************************************************
*** subroutine getqvec
***********************************************************************
SUBROUTINE GETQVEC(L,J,DT,X)
IMPLICIT NONE
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
C--Common block of Pythia
COMMON/PYJETS/N,NPAD,K(23000,5),P(23000,5),V(23000,5)
INTEGER N,NPAD,K
DOUBLE PRECISION P,V
C--time common block
COMMON/TIME/MV(23000,5)
DOUBLE PRECISION MV
C--variables for coherent scattering
COMMON/COHERENT/NSTART,NEND,ALLQS(10000,6),SCATCENTRES(10000,10),
&QSUMVEC(4),QSUM2
INTEGER NSTART,NEND
DOUBLE PRECISION ALLQS,SCATCENTRES,QSUMVEC,QSUM2
C--discard event flag
COMMON/DISC/NDISC,NSTRANGE,NGOOD,errcount,wdisc,DISCARD
LOGICAL DISCARD
INTEGER NDISC,NSTRANGE,NGOOD,errcount
double precision wdisc
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--local variables
INTEGER L,J,COUNTER,COUNTMAX,COUNT2,i
DOUBLE PRECISION XSC,YSC,ZSC,TSC,GETMD,GETTEMP,DT,X,PYR,NEWMOM(4),
&T,PT,MAXT,PHI2,BETA(3),PHI,THETA,GETT,PYP,PI,PT2,GETMS,
&savemom(5),theta2,mb2,pz,kt2,phiq,maxt2,xi,md,shat,pcms2,
&avmom(5)
CHARACTER TYPS
DATA PI/3.141592653589793d0/
DATA COUNTMAX/1000/
IF (J.GT.10000)THEN
discard = .true.
return
ENDIF
COUNTER=0
COUNT2=0
XSC=MV(L,1)+DT*P(L,1)/P(L,4)
YSC=MV(L,2)+DT*P(L,2)/P(L,4)
ZSC=MV(L,3)+DT*P(L,3)/P(L,4)
TSC=MV(L,4)+DT
md = GETMD(XSC,YSC,ZSC,TSC)
call AVSCATCEN(xsc,ysc,zsc,tsc,
&avmom(1),avmom(2),avmom(3),avmom(4),avmom(5))
do 210 i=1,5
savemom(i) = p(l,i)
210 continue
xi = sqrt(max(x**2*p(l,4)**2,p(l,5)**2) - p(l,5)**2)/pyp(l,8)
p(l,1) = xi*p(l,1)
p(l,2) = xi*p(l,2)
p(l,3) = xi*p(l,3)
p(l,4) = max(x*p(l,4),p(l,5))
444 CALL GETSCATTERER(XSC,YSC,ZSC,TSC,
&K(1,2),P(1,1),P(1,2),P(1,3),P(1,4),P(1,5))
MV(1,1)=XSC
MV(1,2)=YSC
MV(1,3)=ZSC
MV(1,4)=TSC
TYPS='Q'
IF(K(1,2).EQ.21)TYPS='G'
shat = avmom(5)**2 + savemom(5)**2 + 2.*(avmom(4)*savemom(4)
& -avmom(1)*savemom(1)-avmom(2)*savemom(2)-avmom(3)*savemom(3))
pcms2 = (shat+savemom(5)**2-avmom(5)**2)**2/(4.*shat)
& -savemom(5)**2
maxt = 4.*pcms2
K(1,1)=13
SCATCENTRES(J,1)=K(1,2)
SCATCENTRES(J,2)=P(1,1)
SCATCENTRES(J,3)=P(1,2)
SCATCENTRES(J,4)=P(1,3)
SCATCENTRES(J,5)=P(1,4)
SCATCENTRES(J,6)=P(1,5)
SCATCENTRES(J,7)=MV(1,1)
SCATCENTRES(J,8)=MV(1,2)
SCATCENTRES(J,9)=MV(1,3)
SCATCENTRES(J,10)=MV(1,4)
C--transform to scattering centre's rest frame and rotate such that parton momentum is in z-direction
BETA(1)=P(1,1)/P(1,4)
BETA(2)=P(1,2)/P(1,4)
BETA(3)=P(1,3)/P(1,4)
CALL PYROBO(L,L,0d0,0d0,-BETA(1),-BETA(2),-BETA(3))
CALL PYROBO(1,1,0d0,0d0,-BETA(1),-BETA(2),-BETA(3))
THETA=PYP(L,13)
PHI=PYP(L,15)
CALL PYROBO(L,L,0d0,-PHI,0d0,0d0,0d0)
CALL PYROBO(1,1,0d0,-PHI,0d0,0d0,0d0)
CALL PYROBO(L,L,-THETA,0d0,0d0,0d0,0d0)
CALL PYROBO(1,1,-THETA,0d0,0d0,0d0,0d0)
C--pick a t from differential scattering cross section
204 T=-GETT(0.d0,MAXT,md)
202 NEWMOM(4)=P(L,4)+T/(2.*p(1,5))
NEWMOM(3)=(T-2.*P(L,5)**2+2.*p(l,4)*NEWMOM(4))/(2.*P(L,3))
PT2=NEWMOM(4)**2-NEWMOM(3)**2-P(L,5)**2
IF(DABS(PT2).LT.1.d-10) PT2=0.d0
IF(T.EQ.0.d0) PT2=0.d0
IF(PT2.LT.0.d0)THEN
T=0.d0
GOTO 202
ENDIF
PT=SQRT(PT2)
PHI2=PYR(0)*2*PI
NEWMOM(1)=PT*COS(PHI2)
NEWMOM(2)=PT*SIN(PHI2)
P(1,1)=NEWMOM(1)-P(L,1)
P(1,2)=NEWMOM(2)-P(L,2)
P(1,3)=NEWMOM(3)-P(L,3)
P(1,4)=NEWMOM(4)-P(L,4)
P(1,5)=0.d0
C--transformation to lab
CALL PYROBO(L,L,THETA,0d0,0d0,0d0,0d0)
CALL PYROBO(1,1,THETA,0d0,0d0,0d0,0d0)
CALL PYROBO(L,L,0d0,PHI,0d0,0d0,0d0)
CALL PYROBO(1,1,0d0,PHI,0d0,0d0,0d0)
CALL PYROBO(L,L,0d0,0d0,BETA(1),BETA(2),BETA(3))
CALL PYROBO(1,1,0d0,0d0,BETA(1),BETA(2),BETA(3))
ALLQS(J,1)=T
ALLQS(J,2)=P(1,1)
ALLQS(J,3)=P(1,2)
ALLQS(J,4)=P(1,3)
ALLQS(J,5)=P(1,4)
QSUMVEC(1)=QSUMVEC(1)+ALLQS(NEND,2)
QSUMVEC(2)=QSUMVEC(2)+ALLQS(NEND,3)
QSUMVEC(3)=QSUMVEC(3)+ALLQS(NEND,4)
QSUMVEC(4)=QSUMVEC(4)+ALLQS(NEND,5)
QSUM2=QSUMVEC(4)**2-QSUMVEC(1)**2-QSUMVEC(2)**2-QSUMVEC(3)**2
IF(QSUM2.GT.0.d0)THEN
QSUMVEC(1)=QSUMVEC(1)-ALLQS(NEND,2)
QSUMVEC(2)=QSUMVEC(2)-ALLQS(NEND,3)
QSUMVEC(3)=QSUMVEC(3)-ALLQS(NEND,4)
QSUMVEC(4)=QSUMVEC(4)-ALLQS(NEND,5)
QSUM2=QSUMVEC(4)**2-QSUMVEC(1)**2-QSUMVEC(2)**2-QSUMVEC(3)**2
IF(COUNTER.GT.COUNTMAX)THEN
write(logfid,*)'GETQVEC unable to find q vector'
ALLQS(J,1)=0.d0
ALLQS(J,2)=0.d0
ALLQS(J,3)=0.d0
ALLQS(J,4)=0.d0
ALLQS(J,5)=0.d0
ELSE
COUNTER=COUNTER+1
GOTO 444
ENDIF
ENDIF
do 211 i=1,5
p(l,i) = savemom(i)
211 continue
END
***********************************************************************
*** subroutine dokinematics
***********************************************************************
SUBROUTINE DOKINEMATICS(L,lold,N1,N2,NEWM,RETRYSPLIT,
& TIME,X,Z,QQBAR)
IMPLICIT NONE
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
C--Common block of Pythia
COMMON/PYJETS/N,NPAD,K(23000,5),P(23000,5),V(23000,5)
INTEGER N,NPAD,K
DOUBLE PRECISION P,V
C--time common block
COMMON/TIME/MV(23000,5)
DOUBLE PRECISION MV
C--factor in front of formation times
COMMON/FTIMEFAC/FTFAC
DOUBLE PRECISION FTFAC
C--colour index common block
COMMON/COLOUR/TRIP(23000),ANTI(23000),COLMAX
INTEGER TRIP,ANTI,COLMAX
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--discard event flag
COMMON/DISC/NDISC,NSTRANGE,NGOOD,errcount,wdisc,DISCARD
LOGICAL DISCARD
INTEGER NDISC,NSTRANGE,NGOOD,errcount
double precision wdisc
C--variables for angular ordering
COMMON/ANGOR/ZA(23000),ZD(23000),THETAA(23000),QQBARD(23000)
DOUBLE PRECISION ZA,ZD,THETAA
LOGICAL QQBARD
C--variables for coherent scattering
COMMON/COHERENT/NSTART,NEND,ALLQS(10000,6),SCATCENTRES(10000,10),
&QSUMVEC(4),QSUM2
INTEGER NSTART,NEND
DOUBLE PRECISION ALLQS,SCATCENTRES,QSUMVEC,QSUM2
C--number of scattering events
COMMON/CHECK/NSCAT,NSCATEFF,NSPLIT,nspliti,nsplitf,nistry,
&nisfail,nfsfail,nfstry,nttot,ntrej
DOUBLE PRECISION NSCAT,NSCATEFF,NSPLIT,nspliti,nsplitf,nistry,
&nisfail,nfsfail,nfstry,nttot,ntrej
C--event weight
COMMON/WEIGHT/EVWEIGHT,sumofweights
double precision EVWEIGHT,sumofweights
C--local variables
INTEGER L,LINE,N1,N2,J,DIR,lold,nold,colmaxold,statold
DOUBLE PRECISION PYR,PI,BETA(3),THETA,PHI,PYP,PHI2,MAXT,T,
&NEWMASS,DELTAM,DM,TTOT,DMLEFT,LAMBDA,TIME,ENDTIME,X,tmp,
&m32,newm2,shat,theta2,z,gettemp,E3new,E4new,p32,p42,p3old,
&newm,mass2,enew,pt2,pt,pl,m12,firsttime,pcms2
CHARACTER*2 TYP
LOGICAL RETRYSPLIT,QQBAR,QQBARDEC,rejectt,redokin,reshuffle
DATA PI/3.141592653589793d0/
if (newm.ne.p(l,5)) then
if (p(l,5).lt.0.d0) then
nistry = nistry+evweight
else
nfstry = nfstry+evweight
endif
endif
IF((N+2*(n2-n1+1)).GT.22990)THEN
write(logfid,*)'event too long for event record'
DISCARD=.TRUE.
RETURN
ENDIF
firsttime = mv(l,5)
redokin = .false.
nttot=nttot+(n2-n1+1)*evweight
newm2=newm
nold=n
colmaxold=colmax
statold=k(l,1)
204 DELTAM=NEWM2-P(L,5)
DMLEFT=DELTAM
TTOT=0.d0
DO 220 J=N1,N2
TTOT=TTOT+ALLQS(J,1)
220 CONTINUE
LINE=L
DO 222 J=N1,N2
C--projectile type
IF(K(LINE,2).EQ.21)THEN
TYP='GC'
IF(PYR(0).LT.0.5)THEN
DIR=1
ELSE
DIR=-1
ENDIF
ELSE
TYP='QQ'
DIR=0
ENDIF
K(1,1)=5
K(1,2)=SCATCENTRES(J,1)
P(1,1)=SCATCENTRES(J,2)
P(1,2)=SCATCENTRES(J,3)
P(1,3)=SCATCENTRES(J,4)
P(1,4)=SCATCENTRES(J,5)
P(1,5)=SCATCENTRES(J,6)
MV(1,1)=SCATCENTRES(J,7)
MV(1,2)=SCATCENTRES(J,8)
MV(1,3)=SCATCENTRES(J,9)
MV(1,4)=SCATCENTRES(J,10)
T=ALLQS(J,1)
if (t.eq.0.d0) then
rejectt = .true.
else
rejectt = .false.
endif
C--transform to c.m.s. and rotate such that parton momentum is in z-direction
BETA(1)=(P(1,1)+p(line,1))/(P(1,4)+p(line,4))
BETA(2)=(P(1,2)+p(line,2))/(P(1,4)+p(line,4))
BETA(3)=(P(1,3)+p(line,3))/(P(1,4)+p(line,4))
IF ((BETA(1).GT.1.d0).OR.(BETA(2).GT.1.d0).OR.(BETA(3).GT.1.d0)
& .or.(sqrt(beta(1)**2+beta(2)**2+beta(3)**2).gt.1.d0))THEN
reshuffle = .false.
else
reshuffle = .true.
endif
205 if (.not.reshuffle) then
BETA(1)=P(1,1)/P(1,4)
BETA(2)=P(1,2)/P(1,4)
BETA(3)=P(1,3)/P(1,4)
CALL PYROBO(LINE,LINE,0d0,0d0,-BETA(1),-BETA(2),-BETA(3))
CALL PYROBO(1,1,0d0,0d0,-BETA(1),-BETA(2),-BETA(3))
THETA=PYP(LINE,13)
PHI=PYP(LINE,15)
CALL PYROBO(LINE,LINE,0d0,-PHI,0d0,0d0,0d0)
CALL PYROBO(1,1,0d0,-PHI,0d0,0d0,0d0)
CALL PYROBO(LINE,LINE,-THETA,0d0,0d0,0d0,0d0)
CALL PYROBO(1,1,-THETA,0d0,0d0,0d0,0d0)
maxt = -2.*p(1,5)*p(line,4)
if (t.lt.maxt) then
t=0.d0
rejectt = .true.
endif
m12 = -p(line,5)**2
203 enew = p(line,4)+t/(2.*p(1,5))
pl = (t+2.*p(line,4)*enew-2.*m12)/(2.*p(line,3))
pt2 = enew**2-pl**2-m12
if (t.eq.0.d0) pt2 = 0.d0
if (dabs(pt2).lt.1.d-8) pt2 = 0.d0
if (pt2.lt.0.d0) then
write(logfid,*)' This should not have happened: pt^2<0!'
write(logfid,*)t,enew,pl,pt2
t = 0.d0
rejectt = .true.
goto 203
endif
pt = sqrt(pt2)
phi2 = pyr(0)*2.*pi
n=n+2
p(n,1)=pt*cos(phi2)
p(n,2)=pt*sin(phi2)
p(n,3)=pl
p(n,4)=enew
p(n,5)=p(line,5)
!---------------------------------
P(N-1,1)=P(1,1)+P(LINE,1)-P(N,1)
P(N-1,2)=P(1,2)+P(LINE,2)-P(N,2)
P(N-1,3)=P(1,3)+P(LINE,3)-P(N,3)
P(N-1,4)=P(1,4)+P(LINE,4)-P(N,4)
mass2 = P(N-1,4)**2-P(N-1,1)**2-P(N-1,2)**2-P(N-1,3)**2
if ((mass2.lt.0.d0).and.(mass2.gt.-1.-6)) mass2=0.d0
if (mass2.lt.0.d0)
& write(logfid,*)'messed up scattering centres mass^2: ',
& mass2,p(1,5)**2
P(N-1,5)=SQRT(mass2)
if (abs(p(n-1,5)-p(1,5)).gt.1.d-6)
& write(logfid,*)'messed up scattering centres mass: ',
& p(n-1,5),p(1,5),p(l,5)
call flush(logfid)
!---------------------------------
! P(N-1,1)=P(1,1)
! P(N-1,2)=P(1,2)
! P(N-1,3)=P(1,3)
! P(N-1,4)=P(1,4)
! P(N-1,5)=P(1,5)
!---------------------------------
else
CALL PYROBO(LINE,LINE,0d0,0d0,-BETA(1),-BETA(2),-BETA(3))
CALL PYROBO(1,1,0d0,0d0,-BETA(1),-BETA(2),-BETA(3))
if ((p(1,4).lt.0.d0).or.(p(line,4).lt.0.d0)) then
CALL PYROBO(1,1,0d0,0d0,BETA(1),BETA(2),BETA(3))
CALL PYROBO(LINE,LINE,0d0,0d0,BETA(1),BETA(2),BETA(3))
reshuffle = .false.
goto 205
endif
THETA=PYP(LINE,13)
PHI=PYP(LINE,15)
CALL PYROBO(LINE,LINE,0d0,-PHI,0d0,0d0,0d0)
CALL PYROBO(1,1,0d0,-PHI,0d0,0d0,0d0)
CALL PYROBO(LINE,LINE,-THETA,0d0,0d0,0d0,0d0)
CALL PYROBO(1,1,-THETA,0d0,0d0,0d0,0d0)
shat = (p(1,4)+p(line,4))**2
p3old = p(line,3)
maxt = -4.*p(line,3)**2
if (t.lt.maxt) then
t=0.d0
rejectt = .true.
ntrej=ntrej+evweight
endif
theta2 = acos(1.d0+t/(2.*p(line,3)**2))
phi2 = pyr(0)*2.*pi
n=n+2
p(n,1)=p(line,3)*sin(theta2)*cos(phi2)
p(n,2)=p(line,3)*sin(theta2)*sin(phi2)
p(n,3)=p(line,3)*cos(theta2)
p(n,4)=p(line,4)
p(n,5)=p(line,5)
!---------------------------------
P(N-1,1)=P(1,1)+P(LINE,1)-P(N,1)
P(N-1,2)=P(1,2)+P(LINE,2)-P(N,2)
P(N-1,3)=P(1,3)+P(LINE,3)-P(N,3)
P(N-1,4)=P(1,4)+P(LINE,4)-P(N,4)
mass2 = P(N-1,4)**2-P(N-1,1)**2-P(N-1,2)**2-P(N-1,3)**2
if ((mass2.lt.0.d0).and.(mass2.gt.-1.-6)) mass2=0.d0
if (mass2.lt.0.d0)
& write(logfid,*)'messed up scattering centres mass^2: ',
& mass2,p(1,5)**2
P(N-1,5)=SQRT(mass2)
if (abs(p(n-1,5)-p(1,5)).gt.1.d-6)
& write(logfid,*)'messed up scattering centres mass: ',
& p(n-1,5),p(1,5),p(l,5)
call flush(logfid)
!---------------------------------
! P(N-1,1)=P(1,1)
! P(N-1,2)=P(1,2)
! P(N-1,3)=P(1,3)
! P(N-1,4)=P(1,4)
! P(N-1,5)=P(1,5)
!---------------------------------
endif
C--outgoing projectile
ZA(N)=1.d0
THETAA(N)=-1.d0
ZD(N)=Z
QQBARD(N)=QQBAR
K(N,1)=K(LINE,1)
K(N,2)=K(LINE,2)
K(N,3)=L
K(N,4)=0
K(N,5)=0
IF(ALLHAD.and.(.not.rejectt))THEN
IF(K(N,2).EQ.21)THEN
IF(DIR.EQ.1)THEN
TRIP(N)=COLMAX+1
ANTI(N)=ANTI(LINE)
ELSE
TRIP(N)=TRIP(LINE)
ANTI(N)=COLMAX+1
ENDIF
ELSEIF(K(N,2).GT.0)THEN
TRIP(N)=COLMAX+1
ANTI(N)=0
ELSE
TRIP(N)=0
ANTI(N)=COLMAX+1
ENDIF
COLMAX=COLMAX+1
ELSE
TRIP(N)=TRIP(LINE)
ANTI(N)=ANTI(LINE)
ENDIF
C--take care of incoming projectile
IF(K(LINE,1).EQ.1)THEN
K(LINE,1)=12
ELSE
K(LINE,1)=14
ENDIF
K(LINE,4)=N-1
K(LINE,5)=N
C--outgoing scattering centre
ZA(N-1)=1.d0
THETAA(N-1)=-1.d0
ZD(N-1)=0.d0
QQBARD(N-1)=.false.
C--temporary status code, will be overwritten later
K(N-1,1)=3
K(N-1,2)=21
K(N-1,3)=0
K(N-1,4)=0
K(N-1,5)=0
IF(ALLHAD.and.(.not.rejectt))THEN
IF((K(N,2).GT.0).AND.(DIR.GE.0))THEN
TRIP(N-1)=TRIP(LINE)
ANTI(N-1)=TRIP(N)
ELSE
TRIP(N-1)=ANTI(N)
ANTI(N-1)=ANTI(LINE)
ENDIF
ELSE
TRIP(N-1)=0
ANTI(N-1)=0
ENDIF
if (reshuffle.and.(dm.gt.0.d0)) then
C--adjust mass and re-shuffle momenta
IF(TTOT.EQ.0.d0)THEN
DM=0.d0
ELSE
if (dmleft.lt.0.d0) then
DM=max(DMLEFT*T/TTOT*1.5d0,dmleft)
else
DM=min(DMLEFT*T/TTOT*1.5d0,dmleft)
endif
ENDIF
TTOT=TTOT-ALLQS(J,1)
newmass = p(n,5)+dm
if (newmass.lt.0.d0) then
m32 = -NEWMASS**2
else
m32 = NEWMASS**2
endif
E3new = (shat + m32 - p(1,5)**2)/(2.d0*sqrt(shat))
E4new = (shat - m32 + p(1,5)**2)/(2.d0*sqrt(shat))
p32 = E3new**2 - m32
p42 = E4new**2 - p(1,5)**2
if ((p32.lt.0.d0).or.(p42.lt.0.d0).or.
& (E3new.lt.0.d0).or.(E4new.lt.0.d0)) then
p32 = 0.d0
p42 = 0.d0
E4new = p(n-1,5)
E3new = sqrt(shat) - E4new
m32 = E3new**2
if ((E3new.lt.0.d0).or.(E4new.lt.0.d0)) then
E3new = p(n,4)
E4new = p(n-1,4)
p32 = p3old**2
p42 = p3old**2
if (p(n,5).lt.0.d0) then
m32 = -p(n,5)**2
else
m32 = p(n,5)**2
endif
endif
endif
p(n,1) = sqrt(p32)*p(n,1)/p3old
p(n,2) = sqrt(p32)*p(n,2)/p3old
p(n,3) = sqrt(p32)*p(n,3)/p3old
p(n,4) = E3new
p(n,5) = sign(sqrt(abs(m32)),newmass)
tmp = p(n,4)**2-p(n,1)**2-p(n,2)**2-p(n,3)**2
if (abs(tmp-m32).gt.1.d-6)
& write(logfid,*) 'Oups, messed up projectiles mass:',
& tmp,m32,p(n,5)
!---------------------------------
p(n-1,1) = sqrt(p42)*p(n-1,1)/p3old
p(n-1,2) = sqrt(p42)*p(n-1,2)/p3old
p(n-1,3) = sqrt(p42)*p(n-1,3)/p3old
p(n-1,4) = E4new
!---------------------------------
! P(N-1,1)=P(1,1)
! P(N-1,2)=P(1,2)
! P(N-1,3)=P(1,3)
! P(N-1,4)=P(1,4)
! P(N-1,5)=P(1,5)
!---------------------------------
tmp = p(n-1,4)**2-p(n-1,1)**2-p(n-1,2)**2-p(n-1,3)**2
& -p(n-1,5)**2
if (abs(tmp).gt.1.d-6)
& write(logfid,*) 'Oups, messed up scattering centres mass:',
& tmp,p3old,p(n-1,1),p(n-1,2),p(n-1,3),p(n-1,4),p(n-1,5)
if ((abs(p(n,1)+p(n-1,1)).gt.1.d-6).or.
& (abs(p(n,2)+p(n-1,2)).gt.1.d-6).or.
& (abs(p(n,3)+p(n-1,3)).gt.1.d-6))
& write(logfid,*) 'Oups, momentum not conserved',
& p(n,1)+p(n-1,1),p(n,2)+p(n-1,2),p(n,3)+p(n-1,3)
endif
C--transformation to lab
CALL PYROBO(N-1,N,THETA,0d0,0d0,0d0,0d0)
CALL PYROBO(LINE,LINE,THETA,0d0,0d0,0d0,0d0)
CALL PYROBO(N-1,N,0d0,PHI,0d0,0d0,0d0)
CALL PYROBO(LINE,LINE,0d0,PHI,0d0,0d0,0d0)
CALL PYROBO(N-1,N,0d0,0d0,BETA(1),BETA(2),BETA(3))
CALL PYROBO(LINE,LINE,0d0,0d0,BETA(1),BETA(2),BETA(3))
CALL PYROBO(1,1,THETA,0d0,0d0,0d0,0d0)
CALL PYROBO(1,1,0d0,PHI,0d0,0d0,0d0)
CALL PYROBO(1,1,0d0,0d0,BETA(1),BETA(2),BETA(3))
if (.not.allhad) then
k(n-1,1)=13
else
IF(SCATRECOIL.AND.(P(N-1,4).GT.1.5*3.*
&GETTEMP(MV(1,1),MV(1,2),MV(1,3),MV(1,4))))THEN
K(N-1,1)=2
ELSE
K(N-1,1)=3
ENDIF
endif
if (rejectt) k(n-1,1)=11
MV(N,4)=MV(1,4)
MV(N-1,4)=MV(1,4)
C--set the production vertices: x_mother + (tprod - tprod_mother) * beta_mother
MV(N-1,1)=MV(line,1)
& +(MV(N-1,4)-MV(line,4))*P(line,1)/max(pyp(line,8),P(line,4))
MV(N-1,2)=MV(line,2)
& +(MV(N-1,4)-MV(line,4))*P(line,2)/max(pyp(line,8),P(line,4))
MV(N-1,3)=MV(line,3)
& +(MV(N-1,4)-MV(line,4))*P(line,3)/max(pyp(line,8),P(line,4))
MV(N, 1)=MV(line,1)
& +(MV(N, 4)-MV(line,4))*P(line,1)/max(pyp(line,8),P(line,4))
MV(N, 2)=MV(line,2)
& +(MV(N, 4)-MV(line,4))*P(line,2)/max(pyp(line,8),P(line,4))
MV(N, 3)=MV(line,3)
& +(MV(N, 4)-MV(line,4))*P(line,3)/max(pyp(line,8),P(line,4))
IF(P(N-1,5).GT.P(1,5))THEN
LAMBDA=1.d0/(FTFAC*0.2*P(N-1,4)/P(N-1,5)**2)
MV(N-1,5)=MV(N-1,4)-LOG(1.d0-PYR(0))/LAMBDA
ELSE
MV(N-1,5)=0.d0
ENDIF
IF(J.LT.N2)THEN
MV(N,5)=SCATCENTRES(J+1,10)
ELSE
IF(P(N,5).GT.0.d0)THEN
IF(DELTAM.EQ.0.d0)THEN
ENDTIME=firsttime
ELSE
IF(X.LT.1.d0)THEN
LAMBDA=1.d0/(FTFAC*P(N,4)*0.2/P(N,5)**2)
ENDTIME=SCATCENTRES(J,10)-LOG(1.d0-PYR(0))/LAMBDA
ELSE
ENDTIME=TIME
ENDIF
ENDIF
MV(N,5)=ENDTIME
ELSE
MV(N,5)=0.d0
ENDIF
ENDIF
MV(LINE,5)=ALLQS(J,6)
if ((.not.redokin).and.(.not.rejectt)) NSCAT=NSCAT+EVWEIGHT
DMLEFT=DMLEFT-(p(n,5)-P(LINE,5))
LINE=N
tmp = abs(p(n,4)**2-p(n,1)**2-p(n,2)**2-p(n,3)**2)-p(n,5)**2
if (abs(tmp).ge.1.d-6)
& write(logfid,*)tmp,j,p(l,5),p(line,5),p(n,5)
222 CONTINUE
if (p(n,5).lt.0.d0) then
nisfail = nisfail+evweight
RETRYSPLIT=.TRUE.
return
endif
if (p(n,5).ne.newm2) then
RETRYSPLIT=.TRUE.
redokin = .true.
nfsfail = nfsfail+evweight
n=nold
colmax=colmaxold
k(l,1)=statold
if (p(l,5).le.0.d0) then
newm2 = 0.d0
else
if (p(l,5).lt.q0) then
if ((newm2.eq.newm).and.(newm.ne.q0+1.d-6)) then
newm2=q0+1.d-6
else
nisfail = nisfail+evweight
RETRYSPLIT=.TRUE.
return
endif
else
newm2=p(l,5)
endif
n2=n1
endif
goto 204
endif
if ((p(n,5).lt.0.d0).or.((p(n,5).gt.0.d0).and.(p(n,5).lt.q0)))
&write(logfid,*)'dokinematics did not reach sensible mass: ',
&p(n,5),newm,p(l,5),newm2
NSCATEFF=NSCATEFF+EVWEIGHT
END
***********************************************************************
*** function getproba
***********************************************************************
DOUBLE PRECISION FUNCTION GETPROBA(QI,QF,QAA,ZAA,EBB,TYPE,
& T1,INS2)
IMPLICIT NONE
C--variables for Sudakov integration
COMMON/SUDAINT/QA,ZA2,EB,T,INSTATE,TYP
DOUBLE PRECISION QA,ZA2,EB,T
CHARACTER*2 TYP
LOGICAL INSTATE
C--local variables
DOUBLE PRECISION QI,QF,QAA,ZAA,EBB,GETSUDAKOV,DERIV,T1
CHARACTER*2 TYPE
LOGICAL INS2
QA=QAA
ZA2=ZAA
EB=EBB
TYP=TYPE
T=T1
INSTATE=INS2
GETPROBA=GETSUDAKOV(QI,QAA,QF,ZAA,EBB,TYPE,T1,INS2)
& *DERIV(QF,1)
END
***********************************************************************
*** function getsudakov
***********************************************************************
DOUBLE PRECISION FUNCTION GETSUDAKOV(QMAX1,QA1,QB1,ZA1,EB1,
& TYPE3,T2,INS)
IMPLICIT NONE
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--variables for Sudakov integration
COMMON/SUDAINT/QA,ZA2,EB,T,INSTATE,TYP
DOUBLE PRECISION QA,ZA2,EB,T
CHARACTER*2 TYP
LOGICAL INSTATE
C--local variables
DOUBLE PRECISION QMAX1,QA1,QB1,ZA1,EB1,TMAX,TB,YSTART,EPSI,
&HFIRST,T2,GETINSUDAFAST,QB2
CHARACTER*2 TYPE3
LOGICAL INS
DATA EPSI/1.d-4/
QB2=QB1
IF(INS)THEN
IF(QB2.LT.Q0) write(logfid,*) 'error: Q < Q0',QB2,QMAX1
IF(QB2.LT.(Q0+1.d-10)) QB2=QB2+1.d-10
ELSE
IF(QB2.LT.Q0) write(logfid,*) 'error: Q < min',QB2,QMAX1
IF(QB2.LT.(Q0+1.d-10)) QB2=QB2+1.d-10
ENDIF
IF(QB2.GE.(QMAX1-1.d-10)) THEN
GETSUDAKOV=1.d0
ELSE
IF(INS)THEN
GETSUDAKOV=GETINSUDAFAST(QB1,QMAX1,TYPE3)
ELSE
QA=QA1
ZA2=ZA1
EB=EB1
TYP=TYPE3
T=T2
INSTATE=.FALSE.
HFIRST=0.01*(QMAX1-QB1)
YSTART=0.d0
CALL ODEINT(YSTART,QB2,QMAX1,EPSI,HFIRST,0.d0,1)
GETSUDAKOV=EXP(-YSTART)
ENDIF
ENDIF
END
***********************************************************************
*** function getinsudakov
***********************************************************************
DOUBLE PRECISION FUNCTION GETINSUDAKOV(QB,QMAX1,TYPE3)
IMPLICIT NONE
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--variables for Sudakov integration
COMMON/SUDAINT/QA,ZA2,EB,T,INSTATE,TYP
DOUBLE PRECISION QA,ZA2,EB,T
CHARACTER*2 TYP
LOGICAL INSTATE
C--local variables
DOUBLE PRECISION QMAX1,QB,QB1,ZA1,EA1,YSTART,EPSI,
&HFIRST
CHARACTER*2 TYPE3
DATA EPSI/1.d-4/
QB1=QB
IF(QB1.LT.Q0) write(logfid,*) 'error: Q < Q0',QB1,QMAX1
IF(QB1.LT.(Q0+1.d-12)) QB1=QB1+1.d-12
IF(QB1.GE.(QMAX1-1.d-12)) THEN
GETINSUDAKOV=1.d0
ELSE
TYP=TYPE3
HFIRST=0.01*(QMAX1-QB1)
YSTART=0.d0
CALL ODEINT(YSTART,QB1,QMAX1,EPSI,HFIRST,0.d0,6)
GETINSUDAKOV=EXP(-YSTART)
ENDIF
END
***********************************************************************
*** function deriv
***********************************************************************
DOUBLE PRECISION FUNCTION DERIV(XVAL,W4)
IMPLICIT NONE
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--variables for splitting function integration
COMMON/INTSPLITF/QQUAD,FM
DOUBLE PRECISION QQUAD,FM
C--variables for Sudakov integration
COMMON/SUDAINT/QA,ZA2,EB,T,INSTATE,TYP
DOUBLE PRECISION QA,ZA2,EB,T
CHARACTER*2 TYP
LOGICAL INSTATE
C--variables for pdf integration
COMMON/PDFINTV/XMAX,Z
DOUBLE PRECISION XMAX,Z
C--variables for cross section integration
COMMON/XSECV/QLOW,MDX
DOUBLE PRECISION QLOW,MDX
C--local variables
INTEGER W4
DOUBLE PRECISION XVAL,GETSPLITI,PI,ALPHAS,GETINSPLITI,
&GETINSUDAFAST,SCATPRIMFUNC,PQQ,PQG,PGG,PGQ,
&MEDDERIV
DATA PI/3.141592653589793d0/
IF(W4.EQ.1)THEN
C--Sudakov integration
IF(INSTATE)THEN
DERIV=2.*GETINSPLITI(XVAL,TYP)/XVAL
ELSE
DERIV=2.*GETSPLITI(QA,XVAL,ZA2,EB,TYP)/XVAL
ENDIF
ELSEIF(W4.EQ.2)THEN
C--P(q->qg) integration
DERIV=(1.+FM)*ALPHAS(XVAL*(1.-XVAL)*QQUAD/1.,LPS)*
& PQQ(XVAL)/(2.*PI)
ELSEIF(W4.EQ.3)THEN
C--P(g->gg) integration
DERIV=(1.+FM)*ALPHAS(XVAL*(1.-XVAL)*QQUAD/1.,LPS)
& *PGG(XVAL)/(2.*PI)
ELSEIF(W4.EQ.4)THEN
C--P(g->qq) integration
DERIV=(1.+FM)*ALPHAS(XVAL*(1-XVAL)*QQUAD/1.,LPS)*
& PQG(XVAL)/(2.*PI)
ELSEIF(W4.EQ.5)THEN
DERIV=EXP(-XVAL)/XVAL
ELSEIF(W4.EQ.6)THEN
DERIV=2.*GETINSPLITI(XVAL,TYP)/XVAL
ELSEIF(W4.EQ.7)THEN
DERIV=2.*GETINSUDAFAST(XVAL,XMAX,'QQ')
& *ALPHAS((1.-Z)*XVAL**2/1.,LPS)
& *PQQ(Z)/(2.*PI*XVAL)
ELSEIF(W4.EQ.8)THEN
DERIV=2.*GETINSUDAFAST(XVAL,XMAX,'GC')
& *ALPHAS((1.-Z)*XVAL**2/1.,LPS)
& *PGQ(Z)/(2.*PI*XVAL)
ELSEIF(W4.EQ.9)THEN
DERIV=2.*GETINSUDAFAST(XVAL,XMAX,'QQ')
& *ALPHAS((1.-Z)*XVAL**2/1.,LPS)
& *PQG(Z)/(2.*PI*XVAL)
ELSEIF(W4.EQ.10)THEN
DERIV=2.*GETINSUDAFAST(XVAL,XMAX,'GC')
& *ALPHAS((1.-Z)*XVAL**2/1.,LPS)*
& *2.*PGG(Z)/(2.*PI*XVAL)
ELSEIF(W4.EQ.11)THEN
DERIV=3.*GETINSPLITI(SCALEFACM*SQRT(XVAL),'GQ')
& *SCATPRIMFUNC(XVAL,MDX)/(2.*XVAL)
ELSEIF(W4.EQ.12)THEN
DERIV=2.*GETINSPLITI(SCALEFACM*SQRT(XVAL),'QG')
& *SCATPRIMFUNC(XVAL,MDX)/(3.*XVAL)
ELSEIF(W4.EQ.13)THEN
DERIV=GETINSUDAFAST(QLOW,SCALEFACM*SQRT(XVAL),'GC')
& *3.*2.*PI*ALPHAS(XVAL+MDX**2,LQCD)**2/(2.*(XVAL+MDX**2)**2)
ELSEIF(W4.EQ.14)THEN
DERIV=GETINSUDAFAST(QLOW,SCALEFACM*SQRT(XVAL),'QQ')
& *2.*2.*PI*ALPHAS(XVAL+MDX**2,LQCD)**2/(3.*(XVAL+MDX**2)**2)
ELSEIF(W4.EQ.21)THEN
DERIV=2.*GETINSUDAFAST(XVAL,XMAX,'QQ')*GETINSPLITI(XVAL,'QQ')
& /XVAL
ELSEIF(W4.EQ.22)THEN
DERIV=2.*GETINSUDAFAST(XVAL,XMAX,'GC')*GETINSPLITI(XVAL,'GQ')
& /XVAL
ELSEIF(W4.EQ.23)THEN
DERIV=2.*GETINSUDAFAST(XVAL,XMAX,'QQ')*GETINSPLITI(XVAL,'QG')
& /XVAL
ELSEIF(W4.EQ.24)THEN
DERIV=2.*GETINSUDAFAST(XVAL,XMAX,'GC')*2.
& *GETINSPLITI(XVAL,'GG')/XVAL
ELSE
DERIV=MEDDERIV(XVAL,W4-100)
ENDIF
END
***********************************************************************
*** function getspliti
***********************************************************************
DOUBLE PRECISION FUNCTION GETSPLITI(QA,QB,ZETA,EB,TYPE1)
IMPLICIT NONE
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--splitting integral
COMMON/SPLITINT/SPLITIGGV(1000,1000),SPLITIQQV(1000,1000),
&SPLITIQGV(1000,1000),QVAL(1000),ZMVAL(1000),QMAX,ZMMIN,NPOINT
INTEGER NPOINT
DOUBLE PRECISION SPLITIGGV,SPLITIQQV,SPLITIQGV,
&QVAL,ZMVAL,QMAX,ZMMIN
C--variables for splitting function integration
COMMON/INTSPLITF/QQUAD,FM
DOUBLE PRECISION QQUAD,FM
C--number of extrapolations in tables
common/extrapolations/ntotspliti,noverspliti,ntotpdf,noverpdf,
&ntotxsec,noverxsec,ntotsuda,noversuda
integer ntotspliti,noverspliti,ntotpdf,noverpdf,
&ntotxsec,noverxsec,ntotsuda,noversuda
C--local variables
INTEGER I,J,LT,QLMAX,ZLMAX,QLINE,ZLINE
DOUBLE PRECISION QA,QB,ZETA,EB,LOW,X1A(2),X2A(2),YA(2,2),Y,
&SPLITINTGG,SPLITINTQG,A,B,YB(2)
CHARACTER*2 TYPE1
ntotspliti=ntotspliti+1
if (qb.gt.qmax) then
noverspliti=noverspliti+1
if (noverspliti.le.25)
& write(logfid,*)'WARNING in getspliti: need to extrapolate: ',
& qb,qmax
endif
C--find boundaries for z integration
IF(ANGORD.AND.(ZETA.NE.1.d0))THEN
LOW=MAX(0.5-0.5*SQRT(1.-Q0**2/QB**2)
& *SQRT(1.-QB**2/EB**2),
& 0.5-0.5*SQRT(1.-4.*QB**2*(1.-ZETA)/(ZETA*QA**2)))
ELSE
LOW=0.5-0.5*SQRT(1.-Q0**2/QB**2)
& *SQRT(1.-QB**2/EB**2)
ENDIF
C--find values in array
QLMAX=INT((QB-QVAL(1))*NPOINT/(QVAL(1000)-QVAL(1))+1)
QLINE=MAX(QLMAX,1)
QLINE=MIN(QLINE,NPOINT)
ZLMAX=INT((LOG(LOW)-LOG(ZMVAL(1)))*NPOINT/
& (LOG(ZMVAL(1000))-LOG(ZMVAL(1)))+1)
ZLINE=MAX(ZLMAX,1)
ZLINE=MIN(ZLINE,NPOINT)
IF((QLINE.GT.999).OR.(ZLINE.GT.999).OR.
& (QLINE.LT.1).OR.(ZLINE.LT.1))THEN
write(logfid,*)'ERROR in GETSPLITI: line number out of bound',
& QLINE,ZLINE
ENDIF
IF((TYPE1.EQ.'GG').OR.(TYPE1.EQ.'GC'))THEN
DO 17 I=1,2
X1A(I)=QVAL(QLINE-1+I)
X2A(I)=ZMVAL(ZLINE-1+I)
DO 16 J=1,2
YA(I,J)=SPLITIGGV(QLINE-1+I,ZLINE-1+J)
16 CONTINUE
17 CONTINUE
DO 30 I=1,2
A=(YA(I,2)-YA(I,1))/(X2A(2)-X2A(1))
B=YA(I,1)-A*X2A(1)
YB(I)=A*LOW+B
30 CONTINUE
IF(X1A(1).EQ.X1A(2))THEN
Y=(YB(1)+YB(2))/2.
ELSE
A=(YB(2)-YB(1))/(X1A(2)-X1A(1))
B=YB(1)-A*X1A(1)
Y=A*QB+B
ENDIF
IF(TYPE1.EQ.'GG')THEN
GETSPLITI=MIN(Y,10.d0)
ELSE
SPLITINTGG=MIN(Y,10.d0)
ENDIF
ENDIF
IF((TYPE1.EQ.'QG').OR.(TYPE1.EQ.'GC'))THEN
DO 19 I=1,2
X1A(I)=QVAL(QLINE-1+I)
X2A(I)=ZMVAL(ZLINE-1+I)
DO 18 J=1,2
YA(I,J)=SPLITIQGV(QLINE-1+I,ZLINE-1+J)
18 CONTINUE
19 CONTINUE
DO 31 I=1,2
A=(YA(I,2)-YA(I,1))/(X2A(2)-X2A(1))
B=YA(I,1)-A*X2A(1)
YB(I)=A*LOW+B
31 CONTINUE
IF(X1A(1).EQ.X1A(2))THEN
Y=(YB(1)+YB(2))/2.
ELSE
A=(YB(2)-YB(1))/(X1A(2)-X1A(1))
B=YB(1)-A*X1A(1)
Y=A*QB+B
ENDIF
IF(TYPE1.EQ.'QG')THEN
GETSPLITI=NF*MIN(Y,10.d0)
ELSE
SPLITINTQG=NF*MIN(Y,10.d0)
ENDIF
ENDIF
IF(TYPE1.EQ.'QQ')THEN
DO 21 I=1,2
X1A(I)=QVAL(QLINE-1+I)
X2A(I)=ZMVAL(ZLINE-1+I)
DO 20 J=1,2
YA(I,J)=SPLITIQQV(QLINE-1+I,ZLINE-1+J)
20 CONTINUE
21 CONTINUE
DO 32 I=1,2
A=(YA(I,2)-YA(I,1))/(X2A(2)-X2A(1))
B=YA(I,1)-A*X2A(1)
YB(I)=A*LOW+B
32 CONTINUE
IF(X1A(1).EQ.X1A(2))THEN
Y=(YB(1)+YB(2))/2.
ELSE
A=(YB(2)-YB(1))/(X1A(2)-X1A(1))
B=YB(1)-A*X1A(1)
Y=A*QB+B
ENDIF
GETSPLITI=MIN(Y,10.d0)
ENDIF
IF(TYPE1.EQ.'GC') GETSPLITI=SPLITINTGG+SPLITINTQG
END
***********************************************************************
*** function getinspliti
***********************************************************************
DOUBLE PRECISION FUNCTION GETINSPLITI(QB,TYPE1)
IMPLICIT NONE
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--local variables
DOUBLE PRECISION QB,LOW,PI,Y,SPLITINTGG,SPLITINTQG,UP,EI
CHARACTER*2 TYPE1
DATA PI/3.141592653589793d0/
C--find boundaries for z integration
UP = 1. - Q0**2/(4.*QB**2)
IF((TYPE1.EQ.'GG').OR.(TYPE1.EQ.'GC'))THEN
LOW=1.d0-UP
IF (UP.LE.LOW) THEN
GETINSPLITI=0.d0
RETURN
ENDIF
Y = 2.* ( LOG(LOG((1.-LOW)*QB**2/LPS**2))
& - LPS**2*EI(LOG((1.-LOW)*QB**2/LPS**2))/QB**2
& + LPS**4*EI(2.*LOG((1.-LOW)*QB**2/LPS**2))/QB**4
& - LPS**6*EI(3.*LOG((1.-LOW)*QB**2/LPS**2))/QB**6
& - LOG(LOG((1.-UP)*QB**2/LPS**2))
& + LPS**2*EI(LOG((1.-UP)*QB**2/LPS**2))/QB**2
& - LPS**4*EI(2.*LOG((1.-UP)*QB**2/LPS**2))/QB**4
& + LPS**6*EI(3.*LOG((1.-UP)*QB**2/LPS**2))/QB**6
& + LOW - LOG(LOW) - UP + LOG(UP) )
& *3.*12.*PI/(2.*PI*(33.-2.*NF))
IF(TYPE1.EQ.'GG')THEN
GETINSPLITI=Y
ELSE
SPLITINTGG=Y
ENDIF
ENDIF
IF((TYPE1.EQ.'QG').OR.(TYPE1.EQ.'GC'))THEN
LOW=0.d0
IF (UP.LE.LOW) THEN
GETINSPLITI=0.d0
RETURN
ENDIF
Y = ( 2.*LPS**6*EI(3.*LOG((1.-LOW)*QB**2/LPS**2))/QB**6
& - 2.*LPS**4*EI(2.*LOG((1.-LOW)*QB**2/LPS**2))/QB**4
& + 2.*LPS**2*EI(LOG((1.-LOW)*QB**2/LPS**2))/QB**2
& - 2.*LPS**6*EI(3.*LOG((1.-UP)*QB**2/LPS**2))/QB**6
& + 2.*LPS**4*EI(2.*LOG((1.-UP)*QB**2/LPS**2))/QB**4
& - 2.*LPS**2*EI(LOG((1.-UP)*QB**2/LPS**2))/QB**2 )
& *12.*PI/(2.*2.*PI*(33.-2.*NF))
IF(TYPE1.EQ.'QG')THEN
GETINSPLITI=NF*Y
ELSE
SPLITINTQG=NF*Y
ENDIF
ENDIF
IF(TYPE1.EQ.'QQ')THEN
LOW=0.d0
IF (UP.LE.LOW) THEN
GETINSPLITI=0.d0
RETURN
ENDIF
Y = ( 2.*LOG(LOG((1.-LOW)*QB**2/LPS**2))
& - 2.*LPS**2*EI(LOG((1.-LOW)*QB**2/LPS**2))/QB**2
& + LPS**4*EI(2.*LOG((1.-LOW)*QB**2/LPS**2))/QB**4
& - 2.*LOG(LOG((1.-UP)*QB**2/LPS**2))
& + 2.*LPS**2*EI(LOG((1.-UP)*QB**2/LPS**2))/QB**2
& - LPS**4*EI(2.*LOG((1.-UP)*QB**2/LPS**2))/QB**4 )
& *4.*12.*PI/(3.*2.*PI*(33.-2.*NF))
GETINSPLITI=Y
ENDIF
IF(TYPE1.EQ.'GQ')THEN
LOW=1.d0-UP
IF (UP.LE.LOW) THEN
GETINSPLITI=0.d0
RETURN
ENDIF
Y = (UP**2/2.-2.*UP+2.*LOG(UP)-LOW**2/2.+2.*LOW- 2.*LOG(LOW))
& *4.*12.*PI/(3.*2.*PI*(33.-2.*NF)*LOG(QB**2/LPS**2))
GETINSPLITI=Y
ENDIF
IF(TYPE1.EQ.'GC') GETINSPLITI=SPLITINTGG+SPLITINTQG
END
***********************************************************************
*** function getpdf
***********************************************************************
DOUBLE PRECISION FUNCTION GETPDF(X,Q,TYP)
IMPLICIT NONE
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
C--pdf common block
COMMON/PDFS/QINQX(2,1000),GINQX(2,1000),QINGX(2,1000),
&GINGX(2,1000)
DOUBLE PRECISION QINQX,GINQX,QINGX,GINGX
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--variables for pdf integration
COMMON/PDFINTV/XMAX,Z
DOUBLE PRECISION XMAX,Z
C--local variables
DOUBLE PRECISION X,Q,QLOW,QHIGH,YSTART,EPSI,HFIRST
CHARACTER*2 TYP
DATA EPSI/1.d-4/
IF((X.LT.0.d0).OR.(X.GT.1.d0).OR.(Q.LT.Q0))THEN
write(logfid,*)'error in GETPDF: parameter out of bound',X,Q
GETPDF=0.d0
RETURN
ENDIF
IF(TYP.EQ.'QQ')THEN
Z=X
XMAX=Q
C--f_q^q
QLOW=MAX(Q0,Q0/(2.*SQRT(1.-X)))
QHIGH=Q
IF((QLOW.GE.QHIGH*(1.d0-1.d-10)).OR.(X.GT.1.d0-1.d-10))THEN
YSTART=0.d0
ELSE
HFIRST=0.01*(QHIGH-QLOW)
YSTART=0.d0
CALL ODEINT(YSTART,QLOW,QHIGH,EPSI,HFIRST,0.d0,7)
ENDIF
GETPDF=YSTART
ELSEIF(TYP.EQ.'GQ')THEN
Z=X
XMAX=Q
C--f_q^g
QLOW=MAX(Q0,MAX(Q0/(2.*SQRT(X)),Q0/(2.*SQRT(1.-X))))
QHIGH=Q
IF((QLOW.GE.QHIGH*(1.d0-1.d-10)).OR.(X.LT.0.d0+1.d-10)
& .OR.(X.GT.1.d0-1.d-10))THEN
YSTART=0.d0
ELSE
HFIRST=0.01*(QHIGH-QLOW)
YSTART=0.d0
CALL ODEINT(YSTART,QLOW,QHIGH,EPSI,HFIRST,0.d0,8)
ENDIF
GETPDF=YSTART
ELSEIF(TYP.EQ.'QG')THEN
Z=X
XMAX=Q
C--f_q^g
QLOW=MAX(Q0,Q0/(2.*SQRT(1.-X)))
QHIGH=Q
IF((QLOW.GE.QHIGH*(1.d0-1.d-10)).OR.(X.GT.1.d0-1.d-10))THEN
YSTART=0.d0
ELSE
HFIRST=0.01*(QHIGH-QLOW)
YSTART=0.d0
CALL ODEINT(YSTART,QLOW,QHIGH,EPSI,HFIRST,0.d0,9)
ENDIF
GETPDF=YSTART
ELSEIF(TYP.EQ.'GG')THEN
Z=X
XMAX=Q
C--f_q^q
QLOW=MAX(Q0,MAX(Q0/(2.*SQRT(X)),Q0/(2.*SQRT(1.-X))))
QHIGH=Q
IF((QLOW.GE.QHIGH*(1.d0-1.d-10)).OR.(X.LT.0.d0+1.d-10)
& .OR.(X.GT.1.d0-1d-10))THEN
YSTART=0.d0
ELSE
HFIRST=0.01*(QHIGH-QLOW)
YSTART=0.d0
CALL ODEINT(YSTART,QLOW,QHIGH,EPSI,HFIRST,0.d0,10)
ENDIF
GETPDF=YSTART
ELSE
write(logfid,*)'error: pdf-type ',TYP,' does not exist'
GETPDF=0.d0
ENDIF
END
***********************************************************************
*** function getpdfxint
***********************************************************************
DOUBLE PRECISION FUNCTION GETPDFXINT(Q,TYP)
IMPLICIT NONE
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
C--pdf common block
COMMON/PDFS/QINQX(2,1000),GINQX(2,1000),QINGX(2,1000),
&GINGX(2,1000)
DOUBLE PRECISION QINQX,GINQX,QINGX,GINGX
C--number of extrapolations in tables
common/extrapolations/ntotspliti,noverspliti,ntotpdf,noverpdf,
&ntotxsec,noverxsec,ntotsuda,noversuda
integer ntotspliti,noverspliti,ntotpdf,noverpdf,
&ntotxsec,noverxsec,ntotsuda,noversuda
C--local variables
INTEGER J,Q2CLOSE,Q2LINE
DOUBLE PRECISION Q,XA(2),YA(2),Y,A,B
CHARACTER*2 TYP
ntotpdf=ntotpdf+1
if (q**2.gt.QINQX(1,1000)) then
noverpdf=noverpdf+1
if (noverpdf.le.25)
& write(logfid,*)'WARNING in getpdfxint: need to extrapolate: ',
& q**2,QINQX(1,1000)
endif
Q2CLOSE=INT((LOG(Q**2)-LOG(QINQX(1,1)))*999.d0/
& (LOG(QINQX(1,1000))-LOG(QINQX(1,1)))+1)
Q2LINE=MAX(Q2CLOSE,1)
Q2LINE=MIN(Q2LINE,999)
IF((Q2LINE.GT.999).OR.(Q2LINE.LT.1))THEN
write(logfid,*)'ERROR in GETPDFXINT: line number out of bound',
& Q2LINE
ENDIF
IF(TYP.EQ.'QQ')THEN
DO 11 J=1,2
XA(J)=QINQX(1,Q2LINE-1+J)
YA(J)=QINQX(2,Q2LINE-1+J)
11 CONTINUE
ELSEIF(TYP.EQ.'GQ')THEN
DO 13 J=1,2
XA(J)=GINQX(1,Q2LINE-1+J)
YA(J)=GINQX(2,Q2LINE-1+J)
13 CONTINUE
ELSEIF(TYP.EQ.'QG')THEN
DO 15 J=1,2
XA(J)=QINGX(1,Q2LINE-1+J)
YA(J)=QINGX(2,Q2LINE-1+J)
15 CONTINUE
ELSEIF(TYP.EQ.'GG')THEN
DO 17 J=1,2
XA(J)=GINGX(1,Q2LINE-1+J)
YA(J)=GINGX(2,Q2LINE-1+J)
17 CONTINUE
ELSE
write(logfid,*)'error in GETPDFXINT: unknown integral type ',TYP
ENDIF
A=(YA(2)-YA(1))/(XA(2)-XA(1))
B=YA(1)-A*XA(1)
Y=A*Q**2+B
GETPDFXINT=Y
END
***********************************************************************
*** subroutine getpdfxintexact
***********************************************************************
DOUBLE PRECISION FUNCTION GETPDFXINTEXACT(Q,TYP)
IMPLICIT NONE
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--variables for pdf integration
COMMON/PDFINTV/XMAX,Z
DOUBLE PRECISION XMAX,Z
C--local variables
DOUBLE PRECISION Q,EPSI,YSTART,HFIRST
CHARACTER*2 TYP
DATA EPSI/1.d-4/
HFIRST=0.01d0
YSTART=0.d0
XMAX=Q
Z=0.d0
IF(TYP.EQ.'QQ')THEN
CALL ODEINT(YSTART,Q0,Q,EPSI,HFIRST,0.d0,21)
ELSEIF(TYP.EQ.'QG')THEN
CALL ODEINT(YSTART,Q0,Q,EPSI,HFIRST,0.d0,23)
ELSEIF(TYP.EQ.'GQ')THEN
CALL ODEINT(YSTART,Q0,Q,EPSI,HFIRST,0.d0,22)
ELSEIF(TYP.EQ.'GG')THEN
CALL ODEINT(YSTART,Q0,Q,EPSI,HFIRST,0.d0,24)
ENDIF
GETPDFXINTEXACT=YSTART
END
***********************************************************************
*** function getxsecint
***********************************************************************
DOUBLE PRECISION FUNCTION GETXSECINT(TM,MD,TYP2)
IMPLICIT NONE
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--cross secttion common block
COMMON/XSECS/INTQ1(1001,101),INTQ2(1001,101),
&INTG1(1001,101),INTG2(1001,101)
DOUBLE PRECISION INTQ1,INTQ2,INTG1,INTG2
C--variables for cross section integration
COMMON/XSECV/QLOW,MDX
DOUBLE PRECISION QLOW,MDX
C--number of extrapolations in tables
common/extrapolations/ntotspliti,noverspliti,ntotpdf,noverpdf,
&ntotxsec,noverxsec,ntotsuda,noversuda
integer ntotspliti,noverspliti,ntotpdf,noverpdf,
&ntotxsec,noverxsec,ntotsuda,noversuda
C--local variables
INTEGER TLINE,TCLOSE,MDCLOSE,MDLINE,I,J
DOUBLE PRECISION TM,X1A(2),X2A(2),YA(2,2),Y,MD,YB(2),A,B
CHARACTER*2 TYP2
ntotxsec=ntotxsec+1
if (tm.gt.intq1(1000,101)) then
noverxsec=noverxsec+1
if (noverpdf.le.25)
& write(logfid,*)'WARNING in getxsecint: need to extrapolate: ',
& tm,intq1(1000,101)
endif
TCLOSE=INT((LOG(TM)-LOG(INTQ1(1,101)))*999.d0/
& (LOG(INTQ1(1000,101))-LOG(INTQ1(1,101)))+1)
TLINE=MAX(TCLOSE,1)
TLINE=MIN(TLINE,999)
MDCLOSE=INT((MD-INTQ1(1001,1))*99.d0/
&(INTQ1(1001,100)-INTQ1(1001,1))+1)
MDLINE=MAX(MDCLOSE,1)
MDLINE=MIN(MDLINE,99)
IF((TLINE.GT.999).OR.(MDLINE.GT.99)
& .OR.(TLINE.LT.1).OR.(MDLINE.LT.1)) THEN
write(logfid,*)'ERROR in GETXSECINT: line number out of bound',
& TLINE,MDLINE
ENDIF
IF(TYP2.EQ.'QA')THEN
C--first quark integral
DO 12 I=1,2
X1A(I)=INTQ1(1001,MDLINE-1+I)
X2A(I)=INTQ1(TLINE-1+I,101)
DO 11 J=1,2
YA(I,J)=INTQ1(TLINE-1+J,MDLINE-1+I)
11 CONTINUE
12 CONTINUE
ELSEIF(TYP2.EQ.'QB')THEN
C--second quark integral
DO 18 I=1,2
X1A(I)=INTQ2(1001,MDLINE-1+I)
X2A(I)=INTQ2(TLINE-1+I,101)
DO 17 J=1,2
YA(I,J)=INTQ2(TLINE-1+J,MDLINE-1+I)
17 CONTINUE
18 CONTINUE
ELSEIF(TYP2.EQ.'GA')THEN
C--first gluon integral
DO 14 I=1,2
X1A(I)=INTG1(1001,MDLINE-1+I)
X2A(I)=INTG1(TLINE-1+I,101)
DO 13 J=1,2
YA(I,J)=INTG1(TLINE-1+J,MDLINE-1+I)
13 CONTINUE
14 CONTINUE
ELSEIF(TYP2.EQ.'GB')THEN
C--second gluon integral
DO 16 I=1,2
X1A(I)=INTG2(1001,MDLINE-1+I)
X2A(I)=INTG2(TLINE-1+I,101)
DO 15 J=1,2
YA(I,J)=INTG2(TLINE-1+J,MDLINE-1+I)
15 CONTINUE
16 CONTINUE
ELSE
write(logfid,*)'error in GETXSECINT: unknown integral type ',
& TYP2
ENDIF
DO 19 I=1,2
A=(YA(I,2)-YA(I,1))/(X2A(2)-X2A(1))
B=YA(I,1)-A*X2A(1)
YB(I)=A*TM+B
19 CONTINUE
IF(X1A(1).EQ.X1A(2))THEN
Y=YB(1)
ELSE
A=(YB(2)-YB(1))/(X1A(2)-X1A(1))
B=YB(1)-A*X1A(1)
Y=A*MD+B
ENDIF
GETXSECINT=Y
END
***********************************************************************
*** function getinsudafast
***********************************************************************
DOUBLE PRECISION FUNCTION GETINSUDAFAST(Q1,Q2,TYP)
IMPLICIT NONE
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--local variables
DOUBLE PRECISION Q1,Q2,GETINSUDARED
CHARACTER*2 TYP
IF(Q2.LE.Q1)THEN
GETINSUDAFAST=1.d0
ELSEIF(Q1.LE.Q0)THEN
GETINSUDAFAST=GETINSUDARED(Q2,TYP)
ELSE
GETINSUDAFAST=GETINSUDARED(Q2,TYP)/GETINSUDARED(Q1,TYP)
ENDIF
IF(GETINSUDAFAST.GT.1.d0) GETINSUDAFAST=1.d0
IF(GETINSUDAFAST.LT.(-1.d-10))THEN
write(logfid,*)'ERROR: GETINSUDAFAST < 0:',
& GETINSUDAFAST,' for',Q1,' ',Q2,' ',TYP
ENDIF
if (getinsudafast.lt.0.d0) getinsudafast = 0.d0
END
***********************************************************************
*** function getinsudared
***********************************************************************
DOUBLE PRECISION FUNCTION GETINSUDARED(Q,TYP2)
IMPLICIT NONE
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--Sudakov common block
COMMON/INSUDA/SUDAQQ(1000,2),SUDAQG(1000,2),SUDAGG(1000,2),
&SUDAGC(1000,2)
DOUBLE PRECISION SUDAQQ,SUDAQG,SUDAGG,SUDAGC
C--number of extrapolations in tables
common/extrapolations/ntotspliti,noverspliti,ntotpdf,noverpdf,
&ntotxsec,noverxsec,ntotsuda,noversuda
integer ntotspliti,noverspliti,ntotpdf,noverpdf,
&ntotxsec,noverxsec,ntotsuda,noversuda
C--local variables
INTEGER QCLOSE,QBIN,I
DOUBLE PRECISION Q,XA(2),YA(2),Y,A,B
CHARACTER*2 TYP2
ntotsuda=ntotsuda+1
if (q.gt.sudaqq(1000,1)) then
noversuda=noversuda+1
if (noversuda.le.25)
& write(logfid,*)'WARNING in getinsudared: need to extrapolate: ',
& q,sudaqq(1000,1)
endif
QCLOSE=INT((LOG(Q)-LOG(SUDAQQ(1,1)))*999.d0
& /(LOG(SUDAQQ(1000,1))-LOG(SUDAQQ(1,1)))+1)
QBIN=MAX(QCLOSE,1)
QBIN=MIN(QBIN,999)
IF((QBIN.GT.999).OR.(QBIN.LT.1)) THEN
write(logfid,*)
& 'ERROR in GETINSUDARED: line number out of bound',QBIN
ENDIF
IF(TYP2.EQ.'QQ')THEN
DO 16 I=1,2
XA(I)=SUDAQQ(QBIN-1+I,1)
YA(I)=SUDAQQ(QBIN-1+I,2)
16 CONTINUE
ELSEIF(TYP2.EQ.'QG')THEN
DO 17 I=1,2
XA(I)=SUDAQG(QBIN-1+I,1)
YA(I)=SUDAQG(QBIN-1+I,2)
17 CONTINUE
ELSEIF(TYP2.EQ.'GG')THEN
DO 18 I=1,2
XA(I)=SUDAGG(QBIN-1+I,1)
YA(I)=SUDAGG(QBIN-1+I,2)
18 CONTINUE
ELSEIF(TYP2.EQ.'GC')THEN
DO 19 I=1,2
XA(I)=SUDAGC(QBIN-1+I,1)
YA(I)=SUDAGC(QBIN-1+I,2)
19 CONTINUE
ELSE
write(logfid,*)'error in GETINSUDARED: unknown type ',TYP2
ENDIF
A=(YA(2)-YA(1))/(XA(2)-XA(1))
B=YA(1)-A*XA(1)
Y=A*Q+B
GETINSUDARED=Y
IF(GETINSUDARED.LT.(-1.d-10))THEN
write(logfid,*) 'ERROR: GETINSUDARED < 0:',GETINSUDARED,Q,TYP2
ENDIF
if (getinsudared.lt.0.d0) getinsudared = 0.d0
END
***********************************************************************
*** function getsscat
***********************************************************************
DOUBLE PRECISION FUNCTION GETSSCAT(EN,px,py,PZ,MP,LW,TYPE1,TYPE2,
& x,y,z,t,mode)
IMPLICIT NONE
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--variables for cross section integration
COMMON/XSECV/QLOW,MDX
DOUBLE PRECISION QLOW,MDX
C--local variables
integer mode
DOUBLE PRECISION UP,EN,LW,SCATPRIMFUNC,CCOL,MP,
&LOW,GETPDFXINT,GETXSECINT,MDEB,pz,pcms2,shat,
&x,y,z,t,getmd,avmom(5),px,py,getmdmin,getmdmax,pproj,psct
CHARACTER TYPE1,TYPE2
IF(TYPE1.EQ.'Q')THEN
CCOL=2./3.
ELSE
CCOL=3./2.
ENDIF
if (mode.eq.0) then
mdeb = getmd(x,y,z,t)
call avscatcen(x,y,z,t,
& avmom(1),avmom(2),avmom(3),avmom(4),avmom(5))
shat = avmom(5)**2 + mp**2 +
& 2.*(avmom(4)*en - avmom(1)*px - avmom(2)*py - avmom(3)*pz)
pcms2 = (shat+mp**2-avmom(5)**2)**2/(4.*shat)-mp**2
up = 4.*pcms2
else
if (mode.eq.1) then
mdeb = getmdmin()
else
mdeb = getmdmax()
endif
call maxscatcen(avmom(1),avmom(2),avmom(3),avmom(4),avmom(5))
psct = sqrt(avmom(1)**2+avmom(2)**2+avmom(3)**2)
pproj = sqrt(px**2+py**2+pz**2)
shat = avmom(5)**2 + mp**2 + 2.*(en*avmom(4) + pproj*psct)
pcms2 = (shat+mp**2-avmom(5)**2)**2/(4.*shat)-mp**2
up = 4.*pcms2
endif
LOW=LW**2
IF(LOW.GT.UP)THEN
GETSSCAT=0.d0
RETURN
ENDIF
IF((TYPE2.EQ.'C').OR.
& ((TYPE1.EQ.'Q').AND.(TYPE2.EQ.'Q')).OR.
& ((TYPE1.EQ.'G').AND.(TYPE2.EQ.'G')))THEN
GETSSCAT=CCOL*(SCATPRIMFUNC(UP,MDEB)-SCATPRIMFUNC(LOW,MDEB))
ELSE
GETSSCAT=0.d0
ENDIF
LOW=Q0**2/SCALEFACM**2
IF(UP.GT.LOW)THEN
IF(TYPE1.EQ.'Q')THEN
IF((TYPE2.EQ.'C').OR.(TYPE2.EQ.'G'))THEN
GETSSCAT=GETSSCAT+GETPDFXINT(SCALEFACM*SQRT(UP),'GQ')
& *3.*SCATPRIMFUNC(UP,MDEB)/2.
GETSSCAT=GETSSCAT-GETXSECINT(UP,MDEB,'QA')
ENDIF
ELSE
IF((TYPE2.EQ.'C').OR.(TYPE2.EQ.'G'))THEN
GETSSCAT=GETSSCAT+CCOL*(SCATPRIMFUNC(UP,MDEB)-
& SCATPRIMFUNC(LOW,MDEB))
& - GETXSECINT(UP,MDEB,'GB')
ENDIF
IF((TYPE2.EQ.'C').OR.(TYPE2.EQ.'Q'))THEN
GETSSCAT=GETSSCAT+2.*GETPDFXINT(SCALEFACM*SQRT(UP),'QG')
& *2.*SCATPRIMFUNC(UP,MDEB)/3.
GETSSCAT=GETSSCAT-2.*GETXSECINT(UP,MDEB,'GA')
ENDIF
ENDIF
ENDIF
IF(GETSSCAT.LT.-1.d-4)
& write(logfid,*) 'error: cross section < 0',GETSSCAT,'for',
& EN,MP,LW,TYPE1,TYPE2,LW**2,UP
GETSSCAT=MAX(GETSSCAT,0.d0)
END
***********************************************************************
*** function getmass
***********************************************************************
DOUBLE PRECISION FUNCTION GETMASS(QBMIN,QBMAX,THETA,EP,flav,
& MAX2,INS,ZDEC,QQBARDEC)
IMPLICIT NONE
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
C--Common block of Pythia
COMMON/PYJETS/N,NPAD,K(23000,5),P(23000,5),V(23000,5)
INTEGER N,NPAD,K
DOUBLE PRECISION P,V
COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
INTEGER MSTU,MSTJ
DOUBLE PRECISION PARU,PARJ
COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
INTEGER MDCY,MDME,KFDP
DOUBLE PRECISION BRAT
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--time common block
COMMON/TIME/MV(23000,5)
DOUBLE PRECISION MV
C--factor in front of alphas argument
COMMON/ALPHASFAC/PTFAC
DOUBLE PRECISION PTFAC
C--local variables
integer flav
DOUBLE PRECISION qbmin,qbmax,theta,ep,max2,zdec,
&q2min,alphmax,alphas,log14,pref,q2max,sudaover,gmin,
&gmax,arg,cand,eps,trueeps,trueval,oest,weight,getinspliti,
- &r,pyr,z,rz,thetanew,r2,pi,pqq,pgg,pqg
+ &r,pyr,z,rz,thetanew,r2,pi,pqq,pgg,pqg,rmin
LOGICAL INS,QQBARDEC
CHARACTER*2 TYPE
DATA PI/3.141592653589793d0/
if(flav.eq.21)then
type='GC'
else
type='QQ'
endif
q2min = q0**2
alphmax = alphas(3.*ptfac*q2min/16.,lps)
log14 = log(0.25)
IF(abs(flav).lt.6)THEN
pref=4.*alphmax/(3.*2.*PI)
ELSE
pref=29.*alphmax/(8.*2.*PI)
ENDIF
-C--check if virtual mass is allowed, return 0.d0 otherwise
+C--check if phase space available, return 0.d0 otherwise
IF((qbmax.LE.QBMIN).OR.(EP.LT.QBMIN)) THEN
getmass=0.d0
ZDEC=0.d0
QQBARDEC=.FALSE.
RETURN
ENDIF
q2max = qbmax**2
- 21 gmax = pref*log(q2min/(4.*q2max))**2
- r=pyr(0)
+! 21 sudaover = exp(-pref*(log(q2min/(4.*q2max))**2 - log14**2))
+! IF(pyr(0).LE.sudaover)THEN
+ 21 if (q2max-qbmin**2.lt.1e-4)then
+ getmass=qbmin
+ zdec=0.5
+ IF(TYPE.EQ.'QQ')THEN
+ QQBARDEC=.FALSE.
+ ELSE
+ IF(PYR(0).LT.PQG(0.5d0)/(PQG(0.5d0)+PGG(0.5d0)))THEN
+ QQBARDEC=.TRUE.
+ ELSE
+ QQBARDEC=.FALSE.
+ ENDIF
+ endif
+ return
+ endif
+ gmax = pref*log(q2min/(4.*q2max))**2
+ if (qbmin.gt.0.d0) then
+ rmin = exp(pref*log(q2min/(4.*qbmin**2))**2-gmax)
+ else
+ rmin = 0.d0
+ endif
+
+ r=pyr(0)*(1.d0-rmin)+rmin
arg=gmax+log(r)
if(arg.lt.0.d0)then
getmass=0.d0
ZDEC=0.d0
QQBARDEC=.FALSE.
RETURN
endif
+! r=pyr(0)
+! gmin = pref*log14**2
+! gmax = pref*log(q2min/(4.*q2max))**2
+! arg = log(r*exp(gmax)+(1.-r)*exp(gmin))
cand = q2min*exp(sqrt(arg/pref))/4.
eps = q2min/(4.*cand)
- if ((cand.lt.q2min**2).or.(cand.lt.qbmin**2)) then
+ if ((cand.lt.q2min).or.(cand.lt.qbmin**2)) then
getmass=0.d0
ZDEC=0.d0
QQBARDEC=.FALSE.
RETURN
endif
IF((CAND.GT.MAX2**2).OR.(CAND.GT.EP**2))THEN
q2max=cand
goto 21
ENDIF
-! if (ins) then
-! trueval=getinspliti(sqrt(cand),type)
-! oest = -2.*pref*log(eps)
- ! weight = trueval/oest
-! else
+ if (ins) then
+ trueval=getinspliti(sqrt(cand),type)
+ oest = -2.*pref*log(eps)
+ weight = trueval/oest
+ else
C--find true z interval
TRUEEPS=0.5-0.5*SQRT(1.-q2min/cand)
& *SQRT(1.-cand/EP**2)
IF(TRUEEPS.LT.EPS)
& WRITE(logfid,*)'error in getmass: true eps < eps',TRUEEPS,EPS
RZ=PYR(0)
z = 1.-eps**rz
if ((z.lt.trueeps).or.(z.gt.(1.-trueeps))) then
weight = 0.
else
if (abs(flav).lt.6)then
- if (ins) then
- trueval = alphas(ptfac*(1.-z)*cand,lps)*pqq(z)/(2.*pi)
- else
trueval = alphas(ptfac*z*(1.-z)*cand,lps)*pqq(z)/(2.*pi)
- endif
oest = 2.*pref/(1.-z)
weight = trueval/oest
else
if (pyr(0).lt.(17./29.)) z = 1.-z
- if (ins)then
- trueval = alphas(ptfac*(1.-z)*cand,lps)
- & *(pgg(z)+pqg(z))/(2.*pi)
- else
trueval = alphas(ptfac*z*(1.-z)*cand,lps)
& *(pgg(z)+pqg(z))/(2.*pi)
- endif
oest = alphmax*(17./(4.*z)+3./(1.-z))/(2.*pi)
weight = trueval/oest
endif
thetanew = sqrt(cand/(z*(1.-z)))/ep
if (angord.and.(theta.gt.0.).and.(thetanew.gt.theta))
& weight = 0.d0
endif
-! endif
+ endif
IF (WEIGHT.GT.1.d0) WRITE(logfid,*)
& 'problem in getmass: weight> 1',
& WEIGHT,TYPE,EPS,TRUEEPS,Z,CAND
R2=PYR(0)
IF(R2.GT.WEIGHT)THEN
q2max=cand
GOTO 21
ELSE
getmass=sqrt(cand)
if (.not.ins) then
ZDEC=Z
IF(abs(flav).lt.6)THEN
QQBARDEC=.FALSE.
ELSE
IF(PYR(0).LT.PQG(Z)/(PQG(Z)+PGG(Z)))THEN
QQBARDEC=.TRUE.
ELSE
QQBARDEC=.FALSE.
ENDIF
ENDIF
endif
ENDIF
END
***********************************************************************
*** function generatez
***********************************************************************
DOUBLE PRECISION FUNCTION GENERATEZ(TI,EA,EPSI,TYPE)
IMPLICIT NONE
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--local variables
DOUBLE PRECISION TI,EA,EPS,PYR,X,R,HELP,R1,EPSI
CHARACTER*2 TYPE
IF(TI.EQ.0.d0)THEN
EPS=EPSI
ELSE
EPS=MAX(0.5-0.5*SQRT(1.-Q0**2/TI)
& *SQRT(1.-TI/EA**2),EPSI)
ENDIF
IF(EPS.GT.0.5)THEN
GENERATEZ=0.5
GOTO 61
ENDIF
60 R=PYR(0)
IF(TYPE.EQ.'QQ')THEN
X=1.-(1.-EPS)*(EPS/(1.-EPS))**R
R=PYR(0)
IF(R.LT.((1.+X**2)/2.))THEN
GENERATEZ=X
ELSE
GOTO 60
ENDIF
ELSEIF(TYPE.EQ.'GG')THEN
X=1./(1.+((1.-EPS)/EPS)**(1.-2.*R))
R=PYR(0)
HELP=((1.-X)/X+X/(1.-X)+X*(1.-X))/(1./(1.-X)+1./X)
IF(R.LT.HELP)THEN
GENERATEZ=X
ELSE
GOTO 60
ENDIF
ELSE
R=PYR(0)*(1.-2.*EPS)+EPS
R1=PYR(0)/2.
HELP=0.5*(R**2+(1.-R)**2)
IF(R1.LT.HELP)THEN
GENERATEZ=R
ELSE
GOTO 60
ENDIF
ENDIF
61 END
***********************************************************************
*** function scatprimfunc
***********************************************************************
DOUBLE PRECISION FUNCTION SCATPRIMFUNC(T,MDEB)
IMPLICIT NONE
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--local variables
DOUBLE PRECISION T,PI,S,EI,ALPHAS,T1,MDEB
DATA PI/3.141592653589793d0/
SCATPRIMFUNC = 2.*PI*(12.*PI)**2*(
& - EI(-LOG((T+MDEB**2)/LQCD**2))/LQCD**2
& - 1./((T+MDEB**2)*LOG((T+MDEB**2)/LQCD**2)))/(33.-2.*NF)**2
END
***********************************************************************
*** function intpqq
***********************************************************************
DOUBLE PRECISION FUNCTION INTPQQ(Z,Q)
IMPLICIT NONE
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--local variables
DOUBLE PRECISION Z,Q
INTPQQ=6.*4.*(-2.*LOG(LOG(Q**2/LPS**2)
& +LOG(1.-Z)))/((33.-2.*NF)*3.)
END
***********************************************************************
*** function intpgglow
***********************************************************************
DOUBLE PRECISION FUNCTION INTPGGLOW(Z,Q)
IMPLICIT NONE
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--local variables
DOUBLE PRECISION Z,Q
INTPGGLOW=6.*3.*(LOG(LOG(Q**2/LPS**2)+LOG(Z)))/(33.-2.*NF)
END
***********************************************************************
*** function intpgghigh
***********************************************************************
DOUBLE PRECISION FUNCTION INTPGGHIGH(Z,Q)
IMPLICIT NONE
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--local variables
DOUBLE PRECISION Z,Q
INTPGGHIGH=-6.*3.*(LOG(LOG(Q**2/LPS**2)+LOG(1.-Z)))/(33.-2.*NF)
END
***********************************************************************
*** function intpqglow
***********************************************************************
DOUBLE PRECISION FUNCTION INTPQGLOW(Z,Q)
IMPLICIT NONE
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--local variables
DOUBLE PRECISION Z,Q,EI
INTPQGLOW=6.*(LPS**2*EI(LOG(Q**2/LPS**2)+LOG(Z))/Q**2
& - 2.*LPS**4*EI(2.*(LOG(Q**2/LPS**2)+LOG(Z)))/Q**4
& + 2.*LPS**6*EI(3.*(LOG(Q**2/LPS**2)+LOG(Z)))/Q**6)/
&((33.-2.*NF)*2.)
END
***********************************************************************
*** function intpqghigh
***********************************************************************
DOUBLE PRECISION FUNCTION INTPQGHIGH(Z,Q)
IMPLICIT NONE
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--local variables
DOUBLE PRECISION Z,Q,EI
INTPQGHIGH=-6.*(LPS**2*EI(LOG(Q**2/LPS**2)+LOG(1.-Z))/Q**2
& - 2.*LPS**4*EI(2.*(LOG(Q**2/LPS**2)+LOG(1.-Z)))/Q**4
& + 2.*LPS**6*EI(3.*(LOG(Q**2/LPS**2)+LOG(1.-Z)))/Q**6)/
&((33.-2.*NF)*2.)
END
***********************************************************************
*** function gett
***********************************************************************
DOUBLE PRECISION FUNCTION GETT(MINT,MAXT,MDEB)
IMPLICIT NONE
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--local variables
DOUBLE PRECISION TMIN,TMAX,MAXI,PYR,R1,R2,ALPHAS,PI,Y,MAXT,
&MDEB,MINT,T
DATA PI/3.141592653589793d0/
TMAX=MAXT+MDEB**2
TMIN=MINT+MDEB**2
IF(TMIN.GT.TMAX) THEN
GETT=0.d0
RETURN
ENDIF
20 R1=PYR(0)
T=TMAX*TMIN/(TMAX+R1*(TMIN-TMAX))
R2=PYR(0)
IF(R2.LT.ALPHAS(T,LQCD)**2/ALPHAS(TMIN,LQCD)**2)THEN
GETT=T-MDEB**2
ELSE
GOTO 20
ENDIF
END
***********************************************************************
*** function ei
***********************************************************************
DOUBLE PRECISION FUNCTION EI(X)
IMPLICIT NONE
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
C--exponential integral for negative arguments
COMMON/EXPINT/EIX(3,1000),VALMAX,NVAL
INTEGER NVAL
DOUBLE PRECISION EIX,VALMAX
C--local variables
INTEGER K,LINE,LMAX
DOUBLE PRECISION X,R,GA,XA(2),YA(2),Y,DY,A,B
DOUBLE PRECISION YSTART,EPSI,HFIRST
DATA EPSI/1.e-5/
IF(DABS(X).GT.VALMAX)
& write(logfid,*)'warning: value out of array in Ei(x)',X,VALMAX
IF(X.GE.0.d0)THEN
LMAX=INT(X*NVAL/VALMAX)
LINE=MAX(LMAX,1)
LINE=MIN(LINE,999)
IF((LINE.GT.999).OR.(LINE.LT.1)) THEN
write(logfid,*)'ERROR in EIX: line number out of bound',LINE
ENDIF
DO 26 K=1,2
XA(K)=EIX(1,LINE-1+K)
YA(K)=EIX(3,LINE-1+K)
26 CONTINUE
A=(YA(2)-YA(1))/(XA(2)-XA(1))
B=YA(1)-A*XA(1)
Y=A*X+B
ELSE
LMAX=INT(-X*NVAL/VALMAX)
LINE=MAX(LMAX,1)
LINE=MIN(LINE,999)
IF((LINE.GT.999).OR.(LINE.LT.1)) THEN
write(logfid,*)'ERROR in EIX: line number out of bound',LINE
ENDIF
DO 27 K=1,2
XA(K)=EIX(1,LINE-1+K)
YA(K)=EIX(2,LINE-1+K)
27 CONTINUE
A=(YA(2)-YA(1))/(XA(2)-XA(1))
B=YA(1)-A*XA(1)
Y=-A*X+B
ENDIF
EI=Y
END
***********************************************************************
*** function pqq
***********************************************************************
DOUBLE PRECISION FUNCTION PQQ(Z)
IMPLICIT NONE
DOUBLE PRECISION Z
PQQ=4.*(1.+Z**2)/(3.*(1.-Z))
END
***********************************************************************
*** function pgq
***********************************************************************
DOUBLE PRECISION FUNCTION PGQ(Z)
IMPLICIT NONE
DOUBLE PRECISION Z
PGQ=4.*(1.+(1.-Z)**2)/(3.*Z)
END
***********************************************************************
*** function pgg
***********************************************************************
DOUBLE PRECISION FUNCTION PGG(Z)
IMPLICIT NONE
DOUBLE PRECISION Z
PGG=3.*((1.-Z)/Z + Z/(1.-Z) + Z*(1.-Z))
END
***********************************************************************
*** function pqg
***********************************************************************
DOUBLE PRECISION FUNCTION PQG(Z)
IMPLICIT NONE
DOUBLE PRECISION Z
PQG=0.5*(Z**2 + (1.-Z)**2)
END
***********************************************************************
*** function alphas
***********************************************************************
DOUBLE PRECISION FUNCTION ALPHAS(T,LAMBDA)
IMPLICIT NONE
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--local variables
DOUBLE PRECISION T,L0,PI,LAMBDA
DATA PI/3.141592653589793d0/
ALPHAS=4.*PI/((11.-2.*NF/3.)*LOG(T/LAMBDA**2))
END
***********************************************************************
*** subroutine splitfncint
***********************************************************************
SUBROUTINE SPLITFNCINT(EMAX)
IMPLICIT NONE
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--splitting integral
COMMON/SPLITINT/SPLITIGGV(1000,1000),SPLITIQQV(1000,1000),
&SPLITIQGV(1000,1000),QVAL(1000),ZMVAL(1000),QMAX,ZMMIN,NPOINT
INTEGER NPOINT
DOUBLE PRECISION SPLITIGGV,SPLITIQQV,SPLITIQGV,
&QVAL,ZMVAL,QMAX,ZMMIN
C--variables for splitting function integration
COMMON/INTSPLITF/QQUAD,FM
DOUBLE PRECISION QQUAD,FM
C--max rapidity
common/rapmax/etamax
double precision etamax
C--local variables
INTEGER NSTEP,I,J
DOUBLE PRECISION EMAX,ZMMAX,EPSI,HFIRST,YSTART,LNZMMIN,
&LNZMMAX,ZM,ZM2,Q,GETMSMAX,avmom(5),shat,pcms2
DATA ZMMAX/0.5/
DATA NSTEP/999/
DATA EPSI/1.d-5/
call maxscatcen(avmom(1),avmom(2),avmom(3),avmom(4),avmom(5))
shat = avmom(5)**2 +
& 2.*emax*(avmom(4)+sqrt(avmom(1)**2+avmom(2)**2+avmom(3)**2))
pcms2 = (shat-avmom(5)**2)**2/(4.*shat)
qmax = sqrt(scalefacm*4.*pcms2)
ZMMIN=Q0/EMAX
LNZMMIN=LOG(ZMMIN)
LNZMMAX=LOG(ZMMAX)
NPOINT=NSTEP
DO 100 I=1,NSTEP+1
Q=(I-1)*(QMAX-Q0)/NSTEP+Q0
QVAL(I)=Q
QQUAD=Q**2
DO 110 J=1,NSTEP+1
ZM=EXP((J-1)*(LNZMMAX-LNZMMIN)/NSTEP+LNZMMIN)
ZMVAL(J)=ZM
IF(Q**2.LT.Q0**2)THEN
ZM2=0.5
ELSE
ZM2=0.5-0.5*SQRT(1.-Q0**2/Q**2)
ENDIF
ZM=MAX(ZM,ZM2)
IF(ZM.EQ.0.5)THEN
SPLITIQQV(I,J)=0.d0
SPLITIGGV(I,J)=0.d0
SPLITIQGV(I,J)=0.d0
ELSE
YSTART=0d0
HFIRST=0.01
FM=0.d0
CALL ODEINT(YSTART,ZM,1.-ZM,EPSI,HFIRST,0d0,2)
SPLITIQQV(I,J)=YSTART
YSTART=0d0
HFIRST=0.01
FM=0.d0
CALL ODEINT(YSTART,ZM,1.-ZM,EPSI,HFIRST,0d0,3)
SPLITIGGV(I,J)=YSTART
YSTART=0d0
HFIRST=0.01
FM=0.d0
CALL ODEINT(YSTART,ZM,1.-ZM,EPSI,HFIRST,0d0,4)
SPLITIQGV(I,J)=YSTART
ENDIF
110 CONTINUE
100 CONTINUE
END
***********************************************************************
*** subroutine pdfint
***********************************************************************
SUBROUTINE PDFINT(EMAX)
IMPLICIT NONE
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--pdf common block
COMMON/PDFS/QINQX(2,1000),GINQX(2,1000),QINGX(2,1000),
&GINGX(2,1000)
DOUBLE PRECISION QINQX,GINQX,QINGX,GINGX
C--variables for pdf integration
COMMON/PDFINTV/XMAX,Z
DOUBLE PRECISION XMAX,Z
C--max rapidity
common/rapmax/etamax
double precision etamax
C--local variables
INTEGER I,J
DOUBLE PRECISION EMAX,Q2,GETPDFXINTEXACT,YSTART,HFIRST,EPSI,
&Q2MAX,DELTAQ2,avmom(5),shat,pcms2
DATA EPSI/1.d-4/
call maxscatcen(avmom(1),avmom(2),avmom(3),avmom(4),avmom(5))
shat = avmom(5)**2 +
& 2.*emax*(avmom(4)+sqrt(avmom(1)**2+avmom(2)**2+avmom(3)**2))
pcms2 = (shat-avmom(5)**2)**2/(4.*shat)
q2max = scalefacm*4.*pcms2
DELTAQ2=LOG(Q2MAX)-LOG(Q0**2)
QINQX(1,1)=Q0**2
GINQX(1,1)=Q0**2
QINGX(1,1)=Q0**2
GINGX(1,1)=Q0**2
QINQX(2,1)=0.d0
GINQX(2,1)=0.d0
QINGX(2,1)=0.d0
GINGX(2,1)=0.d0
DO 12 J=2,1000
Q2 = EXP((J-1)*DELTAQ2/999.d0 + LOG(Q0**2))
QINQX(1,J)=Q2
GINQX(1,J)=Q2
QINGX(1,J)=Q2
GINGX(1,J)=Q2
QINQX(2,J)=GETPDFXINTEXACT(SQRT(Q2),'QQ')
GINQX(2,J)=GETPDFXINTEXACT(SQRT(Q2),'GQ')
QINGX(2,J)=GETPDFXINTEXACT(SQRT(Q2),'QG')
GINGX(2,J)=GETPDFXINTEXACT(SQRT(Q2),'GG')
12 CONTINUE
END
***********************************************************************
*** subroutine xsecint
***********************************************************************
SUBROUTINE XSECINT(EMAX)
IMPLICIT NONE
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--cross secttion common block
COMMON/XSECS/INTQ1(1001,101),INTQ2(1001,101),
&INTG1(1001,101),INTG2(1001,101)
DOUBLE PRECISION INTQ1,INTQ2,INTG1,INTG2
C--variables for cross section integration
COMMON/XSECV/QLOW,MDX
DOUBLE PRECISION QLOW,MDX
C--max rapidity
common/rapmax/etamax
double precision etamax
C--local variables
INTEGER J,K
DOUBLE PRECISION EMAX,TMAX,TMAXMAX,DELTATMAX,YSTART,HFIRST,EPSI,
&GETMSMAX,GETMDMAX,MDMIN,MDMAX,DELTAMD,GETMDMIN,avmom(5),shat,pcms2
DATA EPSI/1.d-4/
call maxscatcen(avmom(1),avmom(2),avmom(3),avmom(4),avmom(5))
shat = avmom(5)**2 +
& 2.*emax*(avmom(4)+sqrt(avmom(1)**2+avmom(2)**2+avmom(3)**2))
pcms2 = (shat-avmom(5)**2)**2/(4.*shat)
tmaxmax = scalefacm*4.*pcms2
DELTATMAX=(LOG(TMAXMAX)-
& LOG(Q0**2*(1.d0+1.d-6)/SCALEFACM**2))/999.d0
MDMIN=GETMDMIN()
MDMAX=MAX(MDMIN,GETMDMAX())
DELTAMD=(MDMAX-MDMIN)/99.d0
DO 12 J=1,1000
TMAX = EXP((J-1)*DELTATMAX
& + LOG(Q0**2*(1.d0+1.d-6)/SCALEFACM**2))
INTQ1(J,101)=TMAX
INTQ2(J,101)=TMAX
INTG1(J,101)=TMAX
INTG2(J,101)=TMAX
DO 13 K=1,100
MDX=MDMIN+(K-1)*DELTAMD
INTQ1(1001,K)=MDX
INTQ2(1001,K)=MDX
INTG1(1001,K)=MDX
INTG2(1001,K)=MDX
IF(TMAX.LT.Q0**2/SCALEFACM**2)THEN
INTQ1(J,K)=0.d0
INTQ2(J,K)=0.d0
INTG1(J,K)=0.d0
INTG2(J,K)=0.d0
ELSE
C--first quark integral
QLOW=Q0
HFIRST=0.01*(TMAX-Q0**2/SCALEFACM**2)
YSTART=0.d0
CALL ODEINT(YSTART,Q0**2/SCALEFACM**2,TMAX,EPSI,HFIRST
& ,0.d0,11)
INTQ1(J,K)=YSTART
C--second quark integral
QLOW=Q0
HFIRST=0.01*(TMAX-Q0**2/SCALEFACM**2)
YSTART=0.d0
CALL ODEINT(YSTART,Q0**2/SCALEFACM**2,TMAX,EPSI,HFIRST
& ,0.d0,14)
INTQ2(J,K)=YSTART
C--first gluon integral
QLOW=Q0
YSTART=0.d0
CALL ODEINT(YSTART,Q0**2/SCALEFACM**2,TMAX,EPSI,HFIRST
& ,0.d0,12)
INTG1(J,K)=YSTART
C--second gluon integral
QLOW=Q0
YSTART=0.d0
CALL ODEINT(YSTART,Q0**2/SCALEFACM**2,TMAX,EPSI,HFIRST
& ,0.d0,13)
INTG2(J,K)=YSTART
ENDIF
13 CONTINUE
12 CONTINUE
END
***********************************************************************
*** function insudaint
***********************************************************************
SUBROUTINE INSUDAINT(EMAX)
IMPLICIT NONE
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--Sudakov common block
COMMON/INSUDA/SUDAQQ(1000,2),SUDAQG(1000,2),SUDAGG(1000,2),
&SUDAGC(1000,2)
DOUBLE PRECISION SUDAQQ,SUDAQG,SUDAGG,SUDAGC
C--max rapidity
common/rapmax/etamax
double precision etamax
C--local variables
INTEGER I
DOUBLE PRECISION QMAX,Q,GETINSUDAKOV,DELTA,EMAX,avmom(5),
&shat,pcms2
call maxscatcen(avmom(1),avmom(2),avmom(3),avmom(4),avmom(5))
shat = avmom(5)**2 +
& 2.*emax*(avmom(4)+sqrt(avmom(1)**2+avmom(2)**2+avmom(3)**2))
pcms2 = (shat-avmom(5)**2)**2/(4.*shat)
qmax = sqrt(scalefacm*4.*pcms2)
DELTA=(LOG(3.*QMAX)-LOG(Q0**2*(1.d0+1.d-6)))/999.d0
DO 22 I=1,1000
Q = EXP((I-1)*DELTA + LOG(Q0**2*(1.d0+1.d-6)))
SUDAQQ(I,1)=Q
SUDAQG(I,1)=Q
SUDAGG(I,1)=Q
SUDAGC(I,1)=Q
SUDAQQ(I,2)=GETINSUDAKOV(Q0,Q,'QQ')
SUDAQG(I,2)=GETINSUDAKOV(Q0,Q,'QG')
SUDAGG(I,2)=GETINSUDAKOV(Q0,Q,'GG')
SUDAGC(I,2)=GETINSUDAKOV(Q0,Q,'GC')
22 CONTINUE
END
***********************************************************************
*** function eixint
***********************************************************************
SUBROUTINE EIXINT
IMPLICIT NONE
C--exponential integral for negative arguments
COMMON/EXPINT/EIX(3,1000),VALMAX,NVAL
INTEGER NVAL
DOUBLE PRECISION EIX,VALMAX
C-local variables
INTEGER I,K
DOUBLE PRECISION X,EPSI,HFIRST,YSTART,EI,GA,R
DATA EPSI/1.d-5/
NVAL=1000
VALMAX=55.
DO 10 I=1,NVAL
X=I*VALMAX/(NVAL*1.d0)
EIX(1,I)=X
C--do negative arguments first
YSTART=0d0
HFIRST=0.01
CALL ODEINT(YSTART,X,1000.d0,EPSI,HFIRST,0.d0,5)
EIX(2,I)=-YSTART
C--now do the positive arguments
IF (X.EQ.0.0) THEN
EI=-1.0D+300
ELSE IF (X.LE.40.0) THEN
EI=1.0D0
R=1.0D0
DO 15 K=1,100
R=R*K*X/(K+1.0D0)**2
EI=EI+R
IF (DABS(R/EI).LE.1.0D-15) GO TO 20
15 CONTINUE
20 GA=0.5772156649015328D0
EI=GA+DLOG(X)+X*EI
ELSE
EI=1.0D0
R=1.0D0
DO 25 K=1,20
R=R*K/X
25 EI=EI+R
EI=DEXP(X)/X*EI
ENDIF
EIX(3,I)=EI
10 CONTINUE
END
***********************************************************************
*** function odeint
***********************************************************************
subroutine odeint(ystart,a,b,eps,h1,hmin,w1)
implicit none
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
C--local variables
integer nmax,nstep,w1
double precision ystart,a,b,eps,h1,hmin,x,h,y,dydx,
&deriv,yscale,hdid,hnew
data nmax/100000/
x = a
y = ystart
h = sign(h1,b-a)
do 20 nstep=1,nmax
dydx = deriv(x,w1)
yscale = abs(y) + abs(h*dydx) + 1.e-25
if (((x + h - b)*h).gt.0.) h = b-x
call rkstepper(x,y,dydx,h,hdid,hnew,yscale,eps,w1)
if ((x - b)*h.ge.0) then
ystart = y
return
endif
h = hnew
if (abs(h).lt.abs(hmin)) then
write(logfid,*)'Error in odeint: stepsize too small',w1
& ,ystart,a,b,h1
return
endif
20 continue
write(logfid,*)'Error in odeint: too many steps',w1
& ,ystart,a,b,h1
end
***********************************************************************
*** function rkstepper
***********************************************************************
subroutine rkstepper(x,y,dydx,htest,hdid,hnew,yscale,eps,w1)
implicit none
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
C--local variables
integer w1
double precision x,y,dydx,htest,hdid,hnew,yscale,eps,
&yhalf,y1,y2,rk4step,dydxhalf,xnew,delta,err,h,safety, powerdown,
&powerup,maxup,maxdown,deriv,fac
logical reject
data powerdown/0.25/
data powerup/0.2/
data safety/0.9/
data maxdown/10./
data maxup/5./
reject = .false.
h = htest
10 xnew = x + h
if (x.eq.xnew) then
write(logfid,*)'Error in rkstepper: step size not significant'
return
endif
yhalf = rk4step(x,y,dydx,h/2.,w1)
dydxhalf = deriv(x+h/2.,w1)
y2 = rk4step(x+h/2.,yhalf,dydxhalf,h/2.,w1)
y1 = rk4step(x,y,dydx,h,w1)
delta = y2-y1
err = abs(delta)/(yscale*eps)
if (err.gt.1.) then
reject = .true.
fac = max(1./maxdown,safety/err**powerdown)
h = h*fac
goto 10
else
if (reject) then
hnew = h
else
fac = min(maxup,safety/err**powerup)
hnew = fac*h
endif
x = xnew
y = y2 + delta/15.
hdid = h
endif
end
***********************************************************************
*** function rk4step
***********************************************************************
double precision function rk4step(x,y,dydx,h,w1)
implicit none
integer w1
double precision x,y,dydx,h,k1,k2,k4,yout,deriv
k1 = h*dydx
k2 = h*deriv(x+h/2.,w1)
k4 = h*deriv(x+h,w1)
yout = y+k1/6.+2.*k2/3.+k4/6.
rk4step = yout
end
***********************************************************************
*** function getdeltat
***********************************************************************
LOGICAL FUNCTION GETDELTAT(LINE,TSTART,DTMAX1,DELTAT)
IMPLICIT NONE
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
C--pythia common block
COMMON/PYJETS/N,NPAD,K(23000,5),P(23000,5),V(23000,5)
INTEGER N,NPAD,K
DOUBLE PRECISION P,V
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--time common block
COMMON/TIME/MV(23000,5)
DOUBLE PRECISION MV
C--max rapidity
common/rapmax/etamax
double precision etamax
C--memory for error message from getdeltat
common/errline/errl
integer errl
C--local variables
INTEGER LINE,I,NNULL
DOUBLE PRECISION DTMAX,SIGMAMAX,NEFFMAX,LINVMAX,PYR,
&R,TOFF,XS,YS,ZS,TS,GETSSCAT,GETMSMAX,GETMDMIN,MSMAX,MDMIN,
&XSTART,YSTART,ZSTART,WEIGHT,MS,MD,NEFF,SIGMA,GETNEFF,
&GETNEFFMAX,GETMS,GETMD,TAU,MDMAX,GETMDMAX,GETNATMDMIN,
&SIGMAMIN,NEFFMIN,TSTART,DTMAX1,DELTAT
CHARACTER PTYPE
LOGICAL STOPNOW
C--initialization
GETDELTAT=.FALSE.
DELTAT=0.D0
DTMAX=DTMAX1
IF(K(LINE,2).EQ.21)THEN
PTYPE='G'
ELSE
PTYPE='Q'
ENDIF
NNULL=0
STOPNOW=.FALSE.
C--check for upper bound from plasma lifetime
IF((TSTART+DTMAX).GE.LTIME)DTMAX=LTIME-TSTART
IF(DTMAX.LT.0.D0) RETURN
C--calculate time relative to production of the considered parton
TOFF=TSTART-MV(LINE,4)
XSTART=MV(LINE,1)+TOFF*P(LINE,1)/P(LINE,4)
YSTART=MV(LINE,2)+TOFF*P(LINE,2)/P(LINE,4)
ZSTART=MV(LINE,3)+TOFF*P(LINE,3)/P(LINE,4)
C--calculate upper limit for density*cross section
SIGMAMAX=GETSSCAT(P(LINE,4),p(line,1),p(line,2),p(line,3),
! & xstart,ystart,-sign(abs(zstart),p(line,3)),zstart+1.d-6)
& P(LINE,5),0.d0,PTYPE,'C',xstart,ystart,zstart,tstart,1)
SIGMAMIN=GETSSCAT(P(LINE,4),p(line,1),p(line,2),p(line,3),
! & xstart,ystart,-sign(abs(zstart),p(line,3)),zstart+1.d-6)
& P(LINE,5),0.d0,PTYPE,'C',xstart,ystart,zstart,tstart,2)
NEFFMAX=GETNEFFMAX()
NEFFMIN=GETNATMDMIN()
LINVMAX=5.d0*MAX(NEFFMIN*SIGMAMAX,NEFFMAX*SIGMAMIN)
DO 333 I=1,1000000
DELTAT=DELTAT-LOG(PYR(0))/LINVMAX
XS=XSTART+DELTAT*P(LINE,1)/P(LINE,4)
YS=YSTART+DELTAT*P(LINE,2)/P(LINE,4)
ZS=ZSTART+DELTAT*P(LINE,3)/P(LINE,4)
TS=TSTART+DELTAT
IF(TS.LT.ZS)THEN
TAU=-1.d0
ELSE
TAU=SQRT(TS**2-ZS**2)
ENDIF
NEFF=GETNEFF(XS,YS,ZS,TS)
IF((TAU.GT.1.d0).AND.(NEFF.EQ.0.d0))THEN
IF(NNULL.GT.4)THEN
STOPNOW=.TRUE.
ELSE
NNULL=NNULL+1
ENDIF
ELSE
NNULL=0
ENDIF
IF((DELTAT.GT.DTMAX).OR.STOPNOW) THEN
DELTAT=DTMAX
RETURN
ENDIF
IF(NEFF.GT.0.d0)THEN
SIGMA=GETSSCAT(P(LINE,4),p(line,1),p(line,2),p(line,3),
& P(LINE,5),0.d0,PTYPE,'C',xs,ys,zs,ts,0)
ELSE
SIGMA=0.d0
ENDIF
WEIGHT=5.d0*NEFF*SIGMA/LINVMAX
IF(WEIGHT.GT.1.d0+1d-6) then
if (line.ne.errl) then
write(logfid,*)'error in GETDELTAT: weight > 1',WEIGHT,
& NEFF*SIGMA/(NEFFMAX*SIGMAMIN),NEFF*SIGMA/(NEFFMIN*SIGMAMAX),
& p(line,4)
errl=line
endif
endif
R=PYR(0)
IF(R.LT.WEIGHT)THEN
GETDELTAT=.TRUE.
RETURN
ENDIF
333 CONTINUE
END
integer function poissonian(lambda)
implicit none
integer n
double precision lambda,disc,p,pyr,u,v,pi
data pi/3.141592653589793d0/
if (lambda.gt.745.d0) then
u = pyr(0);
v = pyr(0);
poissonian =
& int(sqrt(lambda)*sqrt(-2.*log(u))*cos(2.*pi*v)+lambda)
else
disc=exp(-lambda)
p=1.d0
n=0
800 p = p*pyr(0)
if (p.gt.disc) then
n = n+1
goto 800
endif
poissonian=n
endif
end
***********************************************************************
*** function ishadron
***********************************************************************
LOGICAL FUNCTION ISHADRON(ID)
IMPLICIT NONE
C--local variables
INTEGER ID
IF(ABS(ID).LT.100) THEN
ISHADRON=.FALSE.
ELSE
IF(MOD(INT(ABS(ID)/10.),10).EQ.0) THEN
ISHADRON = .FALSE.
ELSE
ISHADRON = .TRUE.
ENDIF
ENDIF
END
***********************************************************************
*** function isdiquark
***********************************************************************
LOGICAL FUNCTION ISDIQUARK(ID)
IMPLICIT NONE
C--local variables
INTEGER ID
IF(ABS(ID).LT.1000) THEN
ISDIQUARK=.FALSE.
ELSE
IF(MOD(INT(ID/10),10).EQ.0) THEN
ISDIQUARK = .TRUE.
ELSE
ISDIQUARK = .FALSE.
ENDIF
ENDIF
END
***********************************************************************
*** function islepton
***********************************************************************
LOGICAL FUNCTION ISLEPTON(ID)
IMPLICIT NONE
C-- local variables
INTEGER ID
IF((ABS(ID).EQ.11).OR.(ABS(ID).EQ.13).OR.(ABS(ID).EQ.15)) THEN
ISLEPTON=.TRUE.
ELSE
ISLEPTON=.FALSE.
ENDIF
END
***********************************************************************
*** function isparton
***********************************************************************
LOGICAL FUNCTION ISPARTON(ID)
IMPLICIT NONE
C--local variables
INTEGER ID
LOGICAL ISDIQUARK
IF((ABS(ID).LT.6).OR.(ID.EQ.21).OR.ISDIQUARK(ID)) THEN
ISPARTON=.TRUE.
ELSE
ISPARTON=.FALSE.
ENDIF
END
***********************************************************************
*** function isprimstring
***********************************************************************
logical function isprimstring(l)
implicit none
COMMON/PYJETS/N,NPAD,K(23000,5),P(23000,5),V(23000,5)
INTEGER N,NPAD,K
DOUBLE PRECISION P,V
C--local variables
integer l
logical isparton
if ((K(l,2).ne.91).and.(K(l,2).ne.92)) then
isprimstring=.false.
return
endif
if ((K(K(l,3),3).eq.0).or.(isparton(K(K(K(l,3),3),2)))) then
isprimstring=.true.
else
isprimstring=.false.
endif
end
***********************************************************************
*** function issecstring
***********************************************************************
logical function issecstring(l)
implicit none
COMMON/PYJETS/N,NPAD,K(23000,5),P(23000,5),V(23000,5)
INTEGER N,NPAD,K
DOUBLE PRECISION P,V
C--local variables
integer l
logical isparton,isprimstring
if ((K(l,2).ne.91).and.(K(l,2).ne.92)) then
issecstring = .false.
return
endif
if (isprimstring(l)) then
issecstring = .false.
return
endif
if (isparton(K(K(K(l,3),3),2))) then
issecstring = .false.
else
issecstring = .true.
endif
end
***********************************************************************
*** function isprimhadron
***********************************************************************
logical function isprimhadron(l)
implicit none
COMMON/PYJETS/N,NPAD,K(23000,5),P(23000,5),V(23000,5)
INTEGER N,NPAD,K
DOUBLE PRECISION P,V
C--local variables
integer l
logical isprimstring,isparton
if (((K(K(l,3),2).EQ.91).OR.(K(K(l,3),2).EQ.92))
& .and.isprimstring(K(l,3))
& .and.(.not.isparton(K(l,2)))) then
isprimhadron=.true.
else
isprimhadron=.false.
endif
if (k(l,1).eq.17) isprimhadron=.true.
end
***********************************************************************
*** function compressevent
***********************************************************************
logical function compressevent(l1)
implicit none
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
COMMON/PYJETS/N,NPAD,K(23000,5),P(23000,5),V(23000,5)
INTEGER N,NPAD,K
DOUBLE PRECISION P,V
C--variables for angular ordering
COMMON/ANGOR/ZA(23000),ZD(23000),THETAA(23000),QQBARD(23000)
DOUBLE PRECISION ZA,ZD,THETAA
LOGICAL QQBARD
C--time common block
COMMON/TIME/MV(23000,5)
DOUBLE PRECISION MV
C--colour index common block
COMMON/COLOUR/TRIP(23000),ANTI(23000),COLMAX
INTEGER TRIP,ANTI,COLMAX
C--local variables
integer l1,i,j,nold,nnew,nstart
nold = n
do 777 i=2,nold
if (((k(i,1).eq.11).or.(k(i,1).eq.12).or.(k(i,1).eq.13)).and.
& (i.ne.l1)) then
nnew = i
goto 778
endif
777 continue
compressevent = .false.
return
778 continue
nstart = nnew
do 779 i=nstart,nold
if (((k(i,1).ne.11).and.(k(i,1).ne.12).and.(k(i,1).ne.13)).or.
& (i.eq.l1)) then
do 780 j=1,5
p(nnew,j)=p(i,j)
v(nnew,j)=v(i,j)
mv(nnew,j)=mv(i,j)
780 continue
trip(nnew)=trip(i)
anti(nnew)=anti(i)
za(nnew)=za(i)
zd(nnew)=zd(i)
thetaa(nnew)=thetaa(i)
qqbard(nnew)=qqbard(i)
k(nnew,1)=k(i,1)
k(nnew,2)=k(i,2)
k(nnew,3)=0
k(nnew,4)=0
k(nnew,5)=0
if (l1.eq.i) l1=nnew
nnew=nnew+1
endif
779 continue
n=nnew-1
if ((nold-n).le.10) then
compressevent = .false.
else
compressevent = .true.
endif
do 781 i=nnew,nold
do 782 j=1,5
k(i,j)=0
p(i,j)=0.d0
v(i,j)=0.d0
mv(i,j)=0.d0
782 continue
trip(i)=0
anti(i)=0
za(i)=0.d0
zd(i)=0.d0
thetaa(i)=0.d0
qqbard(i)=.false.
781 continue
if (n.gt.23000) write(logfid,*)'Error in compressevent: n = ',n
if (l1.gt.n) write(logfid,*)'Error in compressevent: l1 = ',l1
call flush(logfid)
return
end
***********************************************************************
*** subroutine pevrec
***********************************************************************
SUBROUTINE PEVREC(NUM,COL)
C--identifier of file for hepmc output and logfile
implicit none
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
COMMON/PYJETS/N,NPAD,K(23000,5),P(23000,5),V(23000,5)
INTEGER N,NPAD,K
DOUBLE PRECISION P,V
C--variables for angular ordering
COMMON/ANGOR/ZA(23000),ZD(23000),THETAA(23000),QQBARD(23000)
DOUBLE PRECISION ZA,ZD,THETAA
LOGICAL QQBARD
C--time common block
COMMON/TIME/MV(23000,5)
DOUBLE PRECISION MV
C--colour index common block
COMMON/COLOUR/TRIP(23000),ANTI(23000),COLMAX
INTEGER TRIP,ANTI,COLMAX
INTEGER NUM,i
LOGICAL COL
DO 202 I=1,N
V(I,1)=MV(I,1)
V(I,2)=MV(I,2)
V(I,3)=MV(I,3)
V(I,4)=MV(I,4)
V(I,5)=MV(I,5)
IF(COL) write(logfid,*)I,' (',TRIP(I),',',ANTI(I),') [',
&K(I,3),K(I,4),K(I,5),' ] {',K(I,2),K(I,1),' } ',
&ZD(I),THETAA(I)
202 CONTINUE
CALL PYLIST(NUM)
END
***********************************************************************
*** subroutine converttohepmc
***********************************************************************
SUBROUTINE CONVERTTOHEPMC(J,EVNUM,PID,beam1,beam2)
IMPLICIT NONE
COMMON/PYJETS/N,NPAD,K(23000,5),P(23000,5),V(23000,5)
INTEGER N,NPAD,K
DOUBLE PRECISION P,V
COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
INTEGER MSTP,MSTI
DOUBLE PRECISION PARP,PARI
C--Parameter common block
COMMON/PARAM/Q0,LPS,LQCD,LTIME,SCALEFACM,ANGORD,SCATRECOIL,
&ALLHAD,compress,NF
INTEGER NF
DOUBLE PRECISION Q0,LQCD,LTIME,LPS,SCALEFACM
LOGICAL ANGORD,SCATRECOIL,ALLHAD,compress
C--organisation of event record
common/evrecord/nsim,npart,offset,hadrotype,sqrts,collider,hadro,
&shorthepmc,channel,isochannel
integer nsim,npart,offset,hadrotype
double precision sqrts
character*4 collider,channel
character*2 isochannel
logical hadro,shorthepmc
C--production point
common/jetpoint/x0,y0
double precision x0,y0
C--initial pt and virtuality
common/injet/isgluon(2),inpt(2),inmass(2),inphi(2),ineta(2)
integer isgluon
double precision inpt,inmass,inphi,ineta
C--local variables
INTEGER EVNUM,PBARCODE,VBARCODE,CODELIST(25000),I,PID,NSTART,
&NFIRST,NVERTEX,NTOT,J,CODEFIRST
integer intinpt(2),intinmass(2),intinphi(2),intineta
DOUBLE PRECISION mproton,centr,getcentrality,jprodr,phi,pi,
&pdummy,pscatcen,mneutron
LOGICAL ISHADRON,ISDIQUARK,ISPARTON,isprimhadron,isprimstring,
&issecstring
character*2 beam1,beam2
data mproton/0.9383/
data mneutron/0.9396/
DATA PI/3.141592653589793d0/
data pdummy/1.d-6/
5000 FORMAT(A2,I10,I3,3E14.6,2I2,I6,4I2,E14.6)
5100 FORMAT(A2,2E14.6)
! 5200 FORMAT(A2,9I2,4E14.6)
! 5200 FORMAT(A2,2I7,7I2,4E14.6)
5200 FORMAT(A2,6I7,2I2,1I7,4E14.6)
5300 FORMAT(A2,2I2,5E14.6,2I2)
5400 FORMAT(A2,I6,6I2,I6,I2)
5500 FORMAT(A2,I6,I6,5E14.6,3I2,I6,I2)
PBARCODE=0
VBARCODE=0
centr = getcentrality()
jprodr = sqrt(x0**2+y0**2)
if (abs(y0).lt.1.d-8) then
if (x0.gt.0.d0) then
phi = 0.d0
else
phi = pi
endif
else
if (x0.gt.0.d0) then
if (y0.gt.0.d0) then
phi = atan(y0/x0)
else
phi = (3.d0*pi/2.d0) - atan(x0/y0)
endif
else
if (y0.gt.0.d0) then
phi = (pi/2.d0) - atan(x0/y0)
else
phi = pi + atan(y0/x0)
endif
endif
endif
do 140 i=1,2
intinpt(i) = int(inpt(i)*100.)
intinmass(i) = int(inmass(i)*100.)
intinphi(i) = int(inphi(i)*100.)
140 continue
intineta = int(ineta(1)*100.)
if (shorthepmc) then
C--short output
IF(COLLIDER.EQ.'EEJJ')THEN
NVERTEX=3
PBARCODE=5
ELSE
NVERTEX=1
PBARCODE=2
ENDIF
nfirst = 0
do 131 i=1,N
if (((k(i,1).lt.10).or.(k(i,1).eq.17)))
& nfirst = nfirst+1
131 continue
WRITE(J,5000)'E ',EVNUM,-1,0.d0,0.d0,0.d0,0,0,NVERTEX,1,2,0,1,
&PARI(10)
WRITE(J,'(A2,I2,A5)')'N ',1,'"0"'
WRITE(J,'(A)')'U GEV MM'
WRITE(J,5100)'C ',PARI(1)*1.d9,0.d0
WRITE(J,5200)'H ',
& intinpt(1),intinmass(1),intinphi(1),
& intinpt(2),intinmass(2),intinphi(2),
& isgluon(1),isgluon(2),intineta,centr,phi,jprodr,ineta(2)
WRITE(J,5300)'F ',0,0,-1.d0,-1.d0,-1.d0,-1.d0,-1.d0,0,0
C--write out vertex line
IF(COLLIDER.EQ.'EEJJ')THEN
WRITE(J,5400)'V ',-1,0,0,0,0,0,2,1,0
WRITE(J,5500)'P ',1,-11,0.d0,0.d0,sqrts/2.,sqrts/2.,
& 0.00051,2,0,0,-1,0
WRITE(J,5500)'P ',2,11,0.d0,0.d0,-sqrts/2.,sqrts/2.,
& 0.00051,2,0,0,-1,0
WRITE(J,5500)'P ',3,23,0.d0,0.d0,0.d0,sqrts,
& 91.2,2,0,0,-2,0
WRITE(J,5400)'V ',-2,0,0,0,0,0,0,2,0
WRITE(J,5500)'P ',4,PID,sqrts/2.,0.d0,0.d0,sqrts/2.,
& 0.000,2,0,0,-3,0
WRITE(J,5500)'P ',5,-PID,-sqrts/2.,0.d0,0.d0,sqrts/2.,
& 0.000,2,0,0,-3,0
WRITE(J,5400)'V ',-3,0,0,0,0,0,0,NFIRST,0
ELSE
WRITE(J,5400)'V ',-1,0,0,0,0,0,2,NFIRST,0
if (beam1.eq.'p+') then
WRITE(J,5500)'P ',1,2212,0.d0,0.d0,
& sqrt(sqrts**2/4.-mproton**2),sqrts/2.,mproton,2,0,0,-1,0
else
WRITE(J,5500)'P ',1,2212,0.d0,0.d0,
& sqrt(sqrts**2/4.-mneutron**2),sqrts/2.,mneutron,2,0,0,-1,0
endif
if (beam2.eq.'p+') then
WRITE(J,5500)'P ',2,2212,0.d0,0.d0,
& -sqrt(sqrts**2/4.-mproton**2),sqrts/2.,mproton,2,0,0,-1,0
else
WRITE(J,5500)'P ',2,2212,0.d0,0.d0,
& -sqrt(sqrts**2/4.-mneutron**2),sqrts/2.,mneutron,2,0,0,-1,0
endif
ENDIF
C--write out particle lines
do 132 i=1,N
if(((k(i,1).lt.10).or.(k(i,1).eq.17))) then
pbarcode=pbarcode+1
WRITE(J,5500)'P ',PBARCODE,K(I,2),P(I,1),P(I,2),P(I,3),
& P(I,4),P(I,5),1,0,0,0,0
endif
132 continue
else
C--long output
if (hadro) then
C--hadronised events
NFIRST=0
IF(COLLIDER.EQ.'EEJJ')THEN
NVERTEX=3
ELSE
NVERTEX=1
ENDIF
DO 123 I=1,N
IF(K(i,3).ne.0)THEN
NSTART=I
GOTO 124
ENDIF
123 CONTINUE
124 CONTINUE
nstart=0
DO 126 I=NSTART+1,N
IF(isprimhadron(i)) NFIRST=NFIRST+1
IF((ISHADRON(K(I,2)).OR.(ABS(K(I,2)).EQ.15))
& .AND.(K(I,4).NE.0)) NVERTEX=NVERTEX+1
126 CONTINUE
127 CONTINUE
WRITE(J,5000)'E ',EVNUM,-1,0.d0,0.d0,0.d0,0,0,NVERTEX,
&1,2,0,1,PARI(10)
WRITE(J,'(A2,I2,A5)')'N ',1,'"0"'
WRITE(J,'(A)')'U GEV MM'
WRITE(J,5100)'C ',PARI(1)*1.d9,0.d0
WRITE(J,5200)'H ',
& intinpt(1),intinmass(1),intinphi(1),
& intinpt(2),intinmass(2),intinphi(2),
& isgluon(1),isgluon(2),intineta,centr,phi,jprodr,ineta(2)
WRITE(J,5300)'F ',0,0,-1.d0,-1.d0,-1.d0,-1.d0,-1.d0,0,0
C--write out vertex line
IF(COLLIDER.EQ.'EEJJ')THEN
VBARCODE=-3
PBARCODE=5
ELSE
VBARCODE=-1
PBARCODE=2
ENDIF
IF(COLLIDER.EQ.'EEJJ')THEN
WRITE(J,5400)'V ',-1,0,0,0,0,0,2,1,0
WRITE(J,5500)'P ',1,-11,0.d0,0.d0,sqrts/2.,sqrts/2.,
& 0.00051,2,0,0,-1,0
WRITE(J,5500)'P ',2,11,0.d0,0.d0,-sqrts/2.,sqrts/2.,
& 0.00051,2,0,0,-1,0
WRITE(J,5500)'P ',3,23,0.d0,0.d0,0.d0,sqrts,
& 91.2,2,0,0,-2,0
WRITE(J,5400)'V ',-2,0,0,0,0,0,0,2,0
WRITE(J,5500)'P ',4,PID,sqrts/2.,0.d0,0.d0,sqrts/2.,
& 0.000,2,0,0,-3,0
WRITE(J,5500)'P ',5,-PID,-sqrts/2.,0.d0,0.d0,sqrts/2.,
& 0.000,2,0,0,-3,0
WRITE(J,5400)'V ',VBARCODE,0,0,0,0,0,0,NFIRST,0
ELSE
WRITE(J,5400)'V ',-1,0,0,0,0,0,2,NFIRST,0
if (beam1.eq.'p+') then
WRITE(J,5500)'P ',1,2212,0.d0,0.d0,
& sqrt(sqrts**2/4.-mproton**2),sqrts/2.,mproton,2,0,0,-1,0
else
WRITE(J,5500)'P ',1,2212,0.d0,0.d0,
& sqrt(sqrts**2/4.-mneutron**2),sqrts/2.,mneutron,2,0,0,-1,0
endif
if (beam2.eq.'p+') then
WRITE(J,5500)'P ',2,2212,0.d0,0.d0,
& -sqrt(sqrts**2/4.-mproton**2),sqrts/2.,mproton,2,0,0,-1,0
else
WRITE(J,5500)'P ',2,2212,0.d0,0.d0,
& -sqrt(sqrts**2/4.-mneutron**2),sqrts/2.,mneutron,2,0,0,-1,0
endif
ENDIF
CODEFIRST=NFIRST+PBARCODE
C--first write out all particles coming directly from string or cluster decays
DO 125 I=NSTART+1,N
IF(.not.isprimhadron(i))THEN
GOTO 125
ELSE
IF (PBARCODE.EQ.CODEFIRST) GOTO 130
PBARCODE=PBARCODE+1
C--write out particle line
IF(K(I,4).GT.0)THEN
VBARCODE=VBARCODE-1
CODELIST(I)=VBARCODE
WRITE(J,5500)'P ',PBARCODE,K(I,2),P(I,1),P(I,2),P(I,3),
& P(I,4),P(I,5),2,0,0,VBARCODE,0
ELSE
WRITE(J,5500)'P ',PBARCODE,K(I,2),P(I,1),P(I,2),P(I,3),
& P(I,4),P(I,5),1,0,0,0,0
ENDIF
ENDIF
125 CONTINUE
130 CONTINUE
C--now write out all other particles and vertices
DO 129 I=NSTART+1,N
if (isprimhadron(i).or.isprimstring(i)) goto 129
if (isparton(K(i,2))) then
if (ishadron(K(K(i,3),2))) codelist(i)=codelist(K(i,3))
goto 129
endif
if (issecstring(i)) then
codelist(i)=codelist(K(i,3))
goto 129
endif
PBARCODE=PBARCODE+1
IF((K(I,3).NE.K(I-1,3)))THEN
C--write out vertex line
WRITE(J,5400)'V ',CODELIST(K(I,3)),0,0,0,0,0,0,
& K(K(I,3),5)-K(K(I,3),4)+1,0
ENDIF
C--write out particle line
IF(K(I,4).GT.0)THEN
VBARCODE=VBARCODE-1
CODELIST(I)=VBARCODE
WRITE(J,5500)'P ',PBARCODE,K(I,2),P(I,1),P(I,2),P(I,3),
& P(I,4),P(I,5),2,0,0,VBARCODE,0
ELSE
WRITE(J,5500)'P ',PBARCODE,K(I,2),P(I,1),P(I,2),P(I,3),
& P(I,4),P(I,5),1,0,0,0,0
ENDIF
129 CONTINUE
else
C--partonic events
endif
endif
call flush(j)
END
***********************************************************************
*** subroutine printlogo
***********************************************************************
subroutine printlogo(fid)
implicit none
integer fid
write(fid,*)
write(fid,*)' _______________'//
&'__________________________ '
write(fid,*)' | '//
&' | '
write(fid,*)' | JJJJJ EEEEE '//
&' W W EEEEE L | '
write(fid,*)' | J E '//
&' W W E L | '
write(fid,*)' _________________| J EEE '//
&' W W W EEE L |_________________ '
write(fid,*)'| | J J E '//
&' W W W W E L | |'
write(fid,*)'| | JJJ EEEEE '//
&' W W EEEEE LLLLL | |'
write(fid,*)'| |_______________'//
&'__________________________| |'
write(fid,*)'| '//
&' |'
write(fid,*)'| '//
&'this is JEWEL 2.1.0 |'
write(fid,*)'| '//
&' |'
write(fid,*)'| Copyright Korinna C. Zapp (2016)'//
&' [Korinna.Zapp@cern.ch] |'
write(fid,*)'| '//
&' |'
write(fid,*)'| The JEWEL homepage is jewel.hepforge.org '//
&' |'
write(fid,*)'| '//
&' |'
write(fid,*)'| The medium model was partly '//
&'implemented by Jochen Klein |'
write(fid,*)'| [Jochen.Klein@cern.ch]. Raghav '//
&'Kunnawalkam Elayavalli helped with the |'
write(fid,*)'| implementation of the V+jet processes '//
&'[raghav.k.e@cern.ch]. |'
write(fid,*)'| '//
&' |'
write(fid,*)'| Please cite JHEP 1303 (2013) '//
&'080 [arXiv:1212.1599] and optionally |'
write(fid,*)'| EPJC C60 (2009) 617 [arXiv:0804.3568] '//
&'for the physics and arXiv:1311.0048 |'
write(fid,*)'| for the code. The reference for '//
&'V+jet processes is arXiv:xxxx.xxxx. |'
write(fid,*)'| '//
&' |'
write(fid,*)'| JEWEL contains code provided by '//
&'S. Zhang and J. M. Jing |'
write(fid,*)'| (Computation of Special Functions, '//
&'John Wiley & Sons, New York, 1996 and |'
write(fid,*)'| http://jin.ece.illinois.edu) for '//
&'computing the exponential integral Ei(x). |'
write(fid,*)'| '//
&' |'
write(fid,*)'| JEWEL relies heavily on PYTHIA 6'//
&' for the event generation. The modified |'
write(fid,*)'| version of PYTHIA 6.4.25 that is'//
&' shipped with JEWEL is, however, not an |'
write(fid,*)'| official PYTHIA release and must'//
&' not be used for anything else. Please |'
write(fid,*)'| refer to results as "JEWEL+PYTHIA".'//
&' |'
write(fid,*)'| '//
&' |'
write(fid,*)'|_________________________________'//
&'____________________________________________|'
write(fid,*)
write(fid,*)
end
***********************************************************************
*** subroutine printtime
***********************************************************************
subroutine printtime
implicit none
C--identifier of file for hepmc output and logfile
common/hepmcid/hpmcfid,logfid
integer hpmcfid,logfid
C--local variables
integer*4 date(3),time(3)
1000 format (i2.2, '.', i2.2, '.', i4.4, ', ',
& i2.2, ':', i2.2, ':', i2.2 )
call idate(date)
call itime(time)
write(logfid,1000)date,time
end
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Tue, Nov 19, 3:22 PM (1 d, 14 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3804938
Default Alt Text
(198 KB)
Attached To
rJEWELSVN jewelsvn
Event Timeline
Log In to Comment