IAP GITLAB

Commit 4ee587c2 authored by Tanguy Pierog's avatar Tanguy Pierog

upate EPOS 1.99 with latest corrections and add PHOJET, PYTHIA and up-to-date...

upate EPOS 1.99 with latest corrections and add PHOJET, PYTHIA and up-to-date HIJING (problem with gfortran)


git-svn-id: https://devel-ik.fzk.de/svn/mc/epos/branches/root@2454 c7a5e08c-de06-0410-9364-b41cf42a0b17
parent 7c516727
......@@ -44,6 +44,7 @@ epos-omg.f epos-par.f epos-sem.f \
epos-rsh.f epos-qsh.f epos-tim.f \
epos-uti.f epos-xan.f epos-xpr.f \
eposu_no.f \
epos_j.f \
epos-app.f \
models.f \
qgsjet/qgsjet03.f \
......@@ -54,7 +55,14 @@ gheisha/gheisha_2002d.f \
gheisha/gheisha_epos.f \
sibyll/sibyll_21.f \
sibyll/nuclib_21.f \
sibyll/sibyll_epos.f
sibyll/sibyll_epos.f \
phojet/phojet.f \
phojet/phojet_epos.f \
phojet/pythia6115.f \
pythia/pythia_epos.f \
hijing/hijing1.383.f \
hijing/hipyset1.35.f \
hijing/hijing_epos.f
# eposm.f
......@@ -82,7 +90,7 @@ hep:
bin/epos$(SYSTEM): $(OBJS) $(COBJS) $(CXXOBJS)
$(CXX) $(CXXFLAGS) $(OBJS) $(COBJS) $(CXXOBJS) -o $@ $(ROOTLIBS) -l$(G2C)
@(echo "")
@(echo "==> compilation with EPOS, QGSJET and SIBYLL successful!")
@(echo "==> compilation with EPOS, QGSJET, PHOJET, PYTHIA, HIJING, GHEISHA and SIBYLL successful!")
@(echo "==> type bin/epos -h")
@(echo "")
......@@ -116,6 +124,9 @@ dirs:
@if [ ! -d $(LIBDIR)pythia ] ;then \
set -x; mkdir -p $(LIBDIR)pythia; set +x; \
fi
@if [ ! -d $(LIBDIR)phojet ] ;then \
set -x; mkdir -p $(LIBDIR)phojet; set +x; \
fi
@if [ ! -d $(LIBDIR)hijing ] ;then \
set -x; mkdir -p $(LIBDIR)hijing; set +x; \
fi
......
......@@ -2,11 +2,11 @@ c 05.02.2008 Main program and random number generator of epos
c (Do not compile it for CONEX or Corsika)
c-----------------------------------------------------------------------
subroutine aamain
program aamain
c-----------------------------------------------------------------------
include 'epos.inc'
common/photrans/phoele(4),ebeam,noevt
common/photrans/phoele(4),ebeam
save nopeno
call aaset(0)
......@@ -22,7 +22,6 @@ c-----------------------------------------------------------------------
goto 9999
endif
call utpri('aamain',ish,ishini,4)
if(model.ne.1)call IniModel(model)
if(iappl.ne.0)then
if(abs(idprojin).eq.12)then !fake e-A collision with gamma(pi0)-A
......@@ -36,9 +35,7 @@ c-----------------------------------------------------------------------
else
if(ebeam.gt.0.)stop 'ebeam need abs(idproj)=12 for DIS'
endif
noevt=0
do nrebin= 1,abs(noebin)
66 continue
call ainit
! if (nrebin.eq.1) call xini
if(nevent.gt.0)then
......@@ -51,7 +48,7 @@ c-----------------------------------------------------------------------
if(nfull.gt.0)nevent=nfreeze*nfull
if(mod(nevent,nfreeze).ne.0)
& stop'nevent must be a multiple of nfreeze!!!!!!!!! '
if(istore.ne.0) call bstora
if(istore.gt.0) call bstora
do n=1,nevent
if(irewch.eq.1)rewind(ifch)
nfr=mod(n-1,nfreeze)*ispherio
......@@ -60,7 +57,6 @@ c-----------------------------------------------------------------------
if(icotabr.eq.0)call aepos(isign(1,-ninicon)*nin)
if(icocore.ne.0)call IniCon(nin)
enddo
if(noevt.eq.1)goto 66 !fake DIS
77 if(ispherio.ne.0)then
if(nfr.eq.0)write(ifmt,'(a)')
& 'spherio evolution + hadronization ...'
......@@ -75,8 +71,12 @@ c call spherio2(nrevt,nfr)
if(iurqmd.eq.1)call urqmd(n)
call afinal
call xana
if(istore.ge.1.and.istore.le.4) call bstore
if(istore.eq.5) call ustore
if(iappl.ne.3)then
if(istore.ge.1.and.istore.le.3) call bstore
if(istore.eq.4) call lhestore(n)
if(istore.eq.5) call ustore
if(istore.eq.-1) call estore
endif
if(nfr+1.eq.nfreeze.or.ispherio.eq.0)call aseed(1)
enddo
call astati
......@@ -106,7 +106,7 @@ c generates a random number
c-----------------------------------------------------------------------
include 'epos.inc'
double precision dranf
1 rangen=sngl(dranf(dble(rangen)))
1 rangen=sngl(dranf(dble(irandm)))
if(rangen.le.0.)goto 1
if(rangen.ge.1.)goto 1
if(irandm.eq.1)write(ifch,*)'rangen()= ',rangen
......@@ -153,7 +153,7 @@ C=======================================================================
C-----------------------------------------------------------------------
C RAN(DOM NUMBER) GEN(ERATOR) USED IN EPOS
C If calling this function within a DO-loop
C you should use an argument which prevents (dummy) to draw this function
C you should use an argument which prevents (dummy) to draw this function
C outside the loop by an optimizing compiler.
C-----------------------------------------------------------------------
implicit none
......@@ -255,7 +255,7 @@ c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c Convert input seed to EPOS random number seed
c Since input seed and EPOS (from Corsika) seed are different,
c define input seed as : seed=ISEED(3)*1E9+ISEED(2)
c define input seed as : seed=ISEED(3)*1E9+ISEED(2)
c-----------------------------------------------------------------------
IMPLICIT NONE
COMMON /CRRANMA3/CD,CINT,CM,TWOM24,TWOM48,MODCNS
......@@ -290,7 +290,7 @@ c-----------------------------------------------------------------------
common/eporansto2/irndmseq
integer irndmseq
integer iseed(3),iseq,iqq,iseqdum
if(iqq.eq.0)then
irndmseq=iseq
elseif(iqq.eq.-1)then
......@@ -361,7 +361,7 @@ C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
IF ( ISEQ .GT. 0 .AND. ISEQ .LE. KSEQ ) JSEQ = ISEQ
DO IVEC = 1, LENV
UNI = U(I97(JSEQ),JSEQ) - U(J97(JSEQ),JSEQ)
IF ( UNI .LT. 0.D0 ) UNI = UNI + 1.D0
......@@ -430,7 +430,7 @@ C-----------------------------------------------------------------------
SAVE
DATA FIRST / .TRUE. /, IORNDM/11/, JSEQ/1/
C-----------------------------------------------------------------------
IF ( FIRST ) THEN
......@@ -623,7 +623,7 @@ c integer iiseed
c common/eporansto2/irndmseq
c integer irndmseq
c integer iseed(3),iseq,iqq,iseqdum
c
c
c if(iqq.eq.0)then
c call ranfst(seed)
c elseif(iqq.eq.-1)then
......
This diff is collapsed.
......@@ -676,17 +676,20 @@ c ----
mas=0
do l=1,ma
if(la.lt.0)then
ia=iabs(idproj/10)
is=idproj/iabs(idproj)
if(ia.ne.111.and.ia.ne.222.and.ia.ne.333)id=idproj/10*10
if(ia.eq.111.or. ia.eq.222.or. ia.eq.333)id=idproj/10*10+is
if(ia.eq.213)id=1230*is
if(iabs(idproj).lt.20)id=idproj
if(iabs(idproj).lt.20)then
id=idproj
else
ia=iabs(idproj/10)
is=idproj/iabs(idproj)
if(ia.ne.111.and.ia.ne.222.and.ia.ne.333)id=idproj/10*10
if(ia.eq.111.or. ia.eq.222.or. ia.eq.333)id=idproj/10*10+is
if(ia.eq.213)id=1230*is
endif
else
id=1220
if(rangen().le.(la-las)*1./(ma-mas))id=1120
if(id.eq.1120)las=las+1
mas=mas+1
id=1220
if(rangen().le.(la-las)*1./(ma-mas))id=1120
if(id.eq.1120)las=las+1
mas=mas+1
endif
ic1=idtrai(1,id,1)
ic2=idtrai(2,id,1)
......@@ -702,17 +705,20 @@ c ----
mas=0
do l=1,ma
if(la.lt.0)then
ia=iabs(idtarg/10)
is=idtarg/iabs(idtarg)
if(ia.ne.111.and.ia.ne.222.and.ia.ne.333)id=idtarg/10*10
if(ia.eq.111.or. ia.eq.222.or. ia.eq.333)id=idtarg/10*10+is
if(ia.eq.213)id=1230*is
if(iabs(idtarg).lt.20)id=idtarg
if(iabs(idtarg).lt.20)then
id=idtarg
else
ia=iabs(idtarg/10)
is=idtarg/iabs(idtarg)
if(ia.ne.111.and.ia.ne.222.and.ia.ne.333)id=idtarg/10*10
if(ia.eq.111.or. ia.eq.222.or. ia.eq.333)id=idtarg/10*10+is
if(ia.eq.213)id=1230*is
endif
else
id=1220
if(rangen().le.(la-las)*1./(ma-mas))id=1120
if(id.eq.1120)las=las+1
mas=mas+1
id=1220
if(rangen().le.(la-las)*1./(ma-mas))id=1120
if(id.eq.1120)las=las+1
mas=mas+1
endif
ic1=idtrai(1,id,1)
ic2=idtrai(2,id,1)
......@@ -764,7 +770,7 @@ c writes /cptl/
c-----------------------------------------------------------------------
include "epos.inc"
include "epos.incems"
double precision XA(64,3),XB(64,3),BQGS,BMAXQGS,BMAXNEX,BMINNEX
double precision XA(210,3),XB(210,3),BQGS,BMAXQGS,BMAXNEX,BMINNEX
COMMON /Q_QGSNEX1/ XA,XB,BQGS,BMAXQGS,BMAXNEX,BMINNEX
common/nucl3/phi,bimp
common /cncl/xproj(mamx),yproj(mamx),zproj(mamx)
......@@ -773,7 +779,7 @@ c-----------------------------------------------------------------------
double precision bqgs2,bmaxqgs2,bmaxnex2,bminnex2,xan,xbn
common /qgsIInex1/xan(iapmax,3),xbn(iapmax,3)
*,bqgs2,bmaxqgs2,bmaxnex2,bminnex2
common/photrans/phoele(4),ebeam,noevt
common/photrans/phoele(4),ebeam
integer ic(2)
call utpri('conwr ',ish,ishini,6)
......
......@@ -762,7 +762,7 @@ c-----------------------------------------------------------------------
common/dkytab/look(mxlook),cbr(mxdky),mode(5,mxdky)
common/nodcay/nodcay,noeta,nopi0,nonunu,noevol,nohadr
logical nodcay,noeta,nopi0,nonunu,noevol,nohadr
parameter (ndectb=1192)
parameter (ndectb=1193)
real dectab(7,ndectb)
data ((dectab(i,j),i=1,7),j= 1, 18)/
......@@ -1938,7 +1938,7 @@ c *---------------------------------------------
c *-----------k0s()--------------------------
* 20., .68610, 120.,-120., 0., 0., 0.
*, 20.,1.00000, 110., 110., 0., 0., 0./
data ((dectab(i,j),i=1,7),j=1172,1192)/
data ((dectab(i,j),i=1,7),j=1172,ndectb)/
c *-----------k0l-------------------------------
* 320., .2113, 110., 110., 110., 0., 0.
*, 320., .2113, 110., 110., 110., 0., 0.
......@@ -1964,7 +1964,9 @@ c *-----------etac-------------------------------
*, 440., .64 , 220., 110., 110., 0., 0.
*, 440., .76 , 120.,-120., 130.,-130., 0.
*, 440., .88 , 120.,-120., 120.,-120., 0.
*, 440., 1. , 130.,-130., 130.,-130., 0./
*, 440., 1. , 130.,-130., 130.,-130., 0.
c *-----------etac-------------------------------
*,1220., 1. ,1120., 12., -11., 0., 0./
call idresi
......
......@@ -321,6 +321,7 @@ c --- Fix Remnant Excitation
enddo
typevt=0 !ela
if(maproj+matarg.eq.2)then !pp
if(itpr(1).ne.0)then
anintine=anintine+1.
......@@ -329,9 +330,14 @@ c --- Fix Remnant Excitation
anintdiff=anintdiff+1.
if((iep(1).eq.0.and.iet(1).eq.2).or.
& (iet(1).eq.0.and.iep(1).eq.2))anintsdif=anintsdif+1.
if((iep(1).eq.0.and.iet(1).eq.2).or.
& (iet(1).eq.0.and.iep(1).eq.2))typevt=3 !SD
if(iep(1).eq.2.and.iet(1).eq.2)typevt=2 !DD
else
anintine=anintine-1. !diffractive without excitation = elastic
endif
else
typevt=1 !ND
endif
endif
else
......@@ -361,12 +367,22 @@ c --- Fix Remnant Excitation
if(aidif.gt.0.)then
anintdiff=anintdiff+1.
anintine=anintine+1.
if(aidifp.gt.0.5.and.aidift.le.0.5)anintsdif=anintsdif+1.
if(aidifp.gt.0.5.and.aidift.le.0.5)then
anintsdif=anintsdif+1.
typevt=3 !SD
endif
if(aidifp.gt.0.5.and.aidift.gt.0.5)then
typevt=2 !DD
endif
if(ionudi.ne.2)then
if(aidifp.le.0.5.and.aidift.gt.0.5)anintsdif=anintsdif+1.
if(aidifp.le.0.5.and.aidift.gt.0.5)then
anintsdif=anintsdif+1.
typevt=3 !SD
endif
endif
elseif(aiine.gt.0.)then
anintine=anintine+1.
typevt=1 !ND
endif
endif
......@@ -1766,7 +1782,7 @@ c restore x from nuclear splitting
if(npnuc(nuc,1,k).eq.n)then
ipp=irnuc(nuc,1,k)
xpp(ipp)=xpp(ipp)+xxnuc(nuc,1,k)
if(xpp(ipp).ge.1d0)iep(ipp)=0
if(xpp(ipp)-1d0.ge.-1d-10)iep(ipp)=0
xppr(n,k)=xppr(n,k)-xxnuc(nuc,1,k)
xpr(n,k)=xppr(n,k)*xmpr(n,k)
ypr(n,k)=0.5D0*log(xppr(n,k)/xmpr(n,k))
......@@ -1779,7 +1795,7 @@ c restore x from nuclear splitting
if(npnuc(nuc,2,k).eq.n)then
itt=irnuc(nuc,2,k)
xmt(itt)=xmt(itt)+xxnuc(nuc,2,k)
if(xmt(itt).ge.1d0)iet(itt)=0
if(xmt(itt)-1d0.ge.-1d-10)iet(itt)=0
xmpr(n,k)=xmpr(n,k)-xxnuc(nuc,2,k)
xpr(n,k)=xppr(n,k)*xmpr(n,k)
ypr(n,k)=0.5D0*log(xppr(n,k)/xmpr(n,k))
......@@ -6772,6 +6788,7 @@ c limits for momenta
fad=max(1.5,fad)
ptm=p1(5)
amasex=dble(fad*utamnz(jcini,mamod))
fas=2.
if(iremn.eq.3)then
id=idtra(icx,ier,ires,0)
if(ier.eq.0)then
......@@ -6781,7 +6798,6 @@ c limits for momenta
strmas=dble(fas*utamnz(jcfin,mamos))
endif
else
fas=2.
strmas=dble(fas*utamnz(jcfin,mamos))
endif
......
......@@ -1405,7 +1405,12 @@ c pt kink due to split elastic scattering (if zzzz.ne.1.)
write (ifch,*) 'boost of b'
write (ifch,*) 'bx,by,bz,be',bx,by,bz,be
endif
call utrot4(-1,bx,by,bz,ptb(1,ib1+1),ptb(2,ib1+1),ptb(3,ib1+1))
if(abs(bx)+abs(by)+abs(bz).gt.0.)then
call utrot4(-1,bx,by,bz,ptb(1,ib1+1),ptb(2,ib1+1),ptb(3,ib1+1))
else
write(ifmt,*) 'null rot of pt',bx,by,bz
write(ifmt,'(4f12.8)') (ptb(i,ib1+1),i=1,4)
endif
if(ish.ge.8) then
write (ifch,*) 'rot of pt'
write(ifch,'(4f12.8)') (ptb(i,ib1+1),i=1,4)
......@@ -1886,7 +1891,12 @@ c & -pst(2,i)*pst(2,j)-pst(1,i)*pst(1,j)
call utlob4( 1,ax,ay,az,ae,am
& ,ptb(1,ibr),ptb(2,ibr),ptb(3,ibr),ptb(4,ibr))
if(ish.ge.6)write(ifch,'(4f12.8)') (ptb(i,ibr),i=1,4)
call utrot4( 1,bx,by,bz,ptb(1,ibr),ptb(2,ibr),ptb(3,ibr))
if(abs(bx)+abs(by)+abs(bz).gt.0.)then
call utrot4( 1,bx,by,bz,ptb(1,ibr),ptb(2,ibr),ptb(3,ibr))
else
write(ifmt,*) 'null rot of pt (2)',bx,by,bz
write(ifmt,'(4f12.8)') (ptb(i,ibr),i=1,4)
endif
if(ish.ge.6)write(ifch,'(4f12.8)') (ptb(i,ibr),i=1,4)
ax=pst(1,ib)+pst(1,jb)
ay=pst(2,ib)+pst(2,jb)
......@@ -1907,7 +1917,12 @@ c & -pst(2,i)*pst(2,j)-pst(1,i)*pst(1,j)
if(ish.ge.6)write (ifch,*) 'ax,ay,az,ae',ax,ay,az,ae,am
call utlob2(1,ax,ay,az,ae,am,bx,by,bz,be,60)
if(ish.ge.6)write (ifch,*) 'bx,by,bz,be',bx,by,bz,be
call utrot4(-1,bx,by,bz,ptb(1,ibr),ptb(2,ibr),ptb(3,ibr))
if(abs(bx)+abs(by)+abs(bz).gt.0.)then
call utrot4(-1,bx,by,bz,ptb(1,ibr),ptb(2,ibr),ptb(3,ibr))
else
write(ifmt,*) 'null rot of pt (3)',bx,by,bz
write(ifmt,'(4f12.8)') (ptb(i,ibr),i=1,4)
endif
if(ish.ge.6)write(ifch,'(4f12.8)') (ptb(i,ibr),i=1,4)
call utlob4(-1,ax,ay,az,ae,am
& ,ptb(1,ibr),ptb(2,ibr),ptb(3,ibr),ptb(4,ibr))
......
......@@ -305,6 +305,13 @@ c----------------------------------------------------------------------
enddo
p(5)=amass(n)
call utlob2(-1,c(1),c(2),c(3),c(4),c(5),p(1),p(2),p(3),p(4),10)
c protection against very high momentum particle (it can happen very very boosted cluster (which do no really make sense anyway))
if(p(4).ge.0.5*engy)then
nptl=nptlb
iret=1
if(ish.ge.4)write(ifch,*)'Decay skipped (p4 to high) !'
goto 1000
endif
do j=1,5
pptl(j,nptl)=p(j)
enddo
......@@ -342,7 +349,7 @@ c----------------------------------------------------------------------
endif
enddo
if(ish.ge.3)then
if(ish.ge.4)then
write(ifch,*)'decay products:'
call alist('&',nptlb+1,nptl)
if(ish.ge.5)then
......
This diff is collapsed.
......@@ -1950,10 +1950,13 @@ c----------------------------------------------------------------------
xmxms=150d0 !maximum mass for a subcluster
if(ncntpo.eq.1)then
write(ifmt,*)'EPOS used with FUSION option'
if(ish.ge.1)then
print*,'+',('-',i=1,57)
print*,'| ttaus fsgrid nsegce:',ttaus,fsgrid,' ',nsegce
print*,'| delxgr delsgr vocell:',delxgr,delsgr,vocell
print*,'+',('-',i=1,57)
endif
endif
c...compute x,y,z
......
......@@ -766,17 +766,17 @@ c----------------------------------------------------------------------
z0=alpfomi!*epscrw*fscra(ee/egyscr)
if(z0.gt.0.)then
z1=zzp
zpUni=dble(z1**gamfom/z0)*exp(-dble(b*b/delD(1,iclpro,icltar)))
c zpUni=dble(4.*z0*(z1/z0)**1.5)
zzpUni=dble(z1**gamfom/z0)*exp(-dble(b*b/delD(1,iclpro,icltar)))
c zzpUni=dble(4.*z0*(z1/z0)**1.5)
z1=zzt
ztUni=dble(z1**gamfom/z0)*exp(-dble(b*b/delD(1,iclpro,icltar)))
c ztUni=dble(4.*z0*(z1/z0)**1.5)
zztUni=dble(z1**gamfom/z0)*exp(-dble(b*b/delD(1,iclpro,icltar)))
c zztUni=dble(4.*z0*(z1/z0)**1.5)
else
zpUni=0d0
ztUni=0d0
zzpUni=0d0
zztUni=0d0
endif
if(ish.ge.6)write(ifch,*)' GfomPar :',zpUni,ztUni
if(ish.ge.6)write(ifch,*)' GfomPar :',zzpUni,zztUni
call utprjx('GfomPar ',ish,ishini,6)
end
......@@ -2995,7 +2995,7 @@ c if(i.ge.2)imax(i)=imax(i)*2
if(b.ge.1.5)imax(i)=2 !max(2,imax(i)/2)
imax(i)=min(30,imax(i))
if(i.gt.imax0)then
if((zpUni*ztUni.lt.1.d-6)
if((zzpUni*zztUni.lt.1.d-6)
& .or.(xp.lt.0.1d0.and.xm.lt.0.1d0))then
imax(i)=0
else
......@@ -3121,7 +3121,7 @@ c if(i.ge.2)imax(i)=imax(i)*2
if(b.ge.1.5)imax(i)=max(2,imax(i)/2)
imax(i)=min(30,imax(i))
if(i.gt.imax0)then
if((zpUni*ztUni.lt.1.d-6)
if((zzpUni*zztUni.lt.1.d-6)
& .or.(xp.lt.0.1d0.and.xm.lt.0.1d0))then
imax(i)=0
else
......@@ -3444,9 +3444,9 @@ c common /ar3/ x1(7),a1(7)
GbetpUni(imax0+1+i)=utgam2(HbetpUni(imax0+1+i))
GbetpUni(imax0+3+i)=utgam2(HbetpUni(imax0+3+i))
GbetpUni(imax0+5+i)=utgam2(HbetpUni(imax0+5+i))
HalpUni(imax0+1+i)=ztUni*alpUni(i,1)
HalpUni(imax0+3+i)=zpUni*alpUni(i,1)
HalpUni(imax0+5+i)=zpUni*ztUni*alpUni(i,1)
HalpUni(imax0+1+i)=zztUni*alpUni(i,1)
HalpUni(imax0+3+i)=zzpUni*alpUni(i,1)
HalpUni(imax0+5+i)=zzpUni*zztUni*alpUni(i,1)
enddo
w=0.d0
......@@ -3978,7 +3978,7 @@ c Calculation by numeric integration :
Df=Df+dble(alp)*xp**dble(bet)*xm**dble(betp)
enddo
call GfomPar(b,s)
G2=Df*(1.d0+ztUni*xp**betfom)*(1.d0+zpUni*xm**betfom)
G2=Df*(1.d0+zztUni*xp**betfom)*(1.d0+zzpUni*xm**betfom)
w=w+dble(a1(j))*PoInU(xp,xm,s,b)*G2
enddo
enddo
......@@ -4048,7 +4048,7 @@ c Calculation by numeric integration :
Df=Df+alp*xp**dble(bet)*xm**dble(betp)
enddo
call GfomPar(b,s)
G2=Df*(1.d0+ztUni*xp**betfom)*(1.d0+zpUni*xm**betfom)
G2=Df*(1.d0+zztUni*xp**betfom)*(1.d0+zzpUni*xm**betfom)
w=w+dble(a1(j))*PoInU(xp,xm,s,b)*G2
enddo
enddo
......@@ -4267,7 +4267,7 @@ c Calculation by numeric integration :
Df=Df+alp*xp**dble(bet)*xm**dble(betp)
enddo
call GfomPar(b,s)
Df=Df*(1.d0+ztUni*xp**betfom)*(1.d0+zpUni*xm**betfom)
Df=Df*(1.d0+zztUni*xp**betfom)*(1.d0+zzpUni*xm**betfom)
w=w+dble(a1(j))*Df*PoInU(xp,xm,s,b)
enddo
enddo
......@@ -4401,7 +4401,7 @@ c----------------------------------------------------------------------
yp=0.d0
if(xm.ne.0.d0)yp=0.5d0*log(xp/xm)
PomIncUnit=om1(xh,yp,b)*PoInU(xp,xm,engy*engy,b)
& *(1.d0+ztUni*xp**betfom)*(1.d0+zpUni*xm**betfom)
& *(1.d0+zztUni*xp**betfom)*(1.d0+zzpUni*xm**betfom)
return
end
......
......@@ -45,7 +45,7 @@ C
C**********************************************************************
include 'epos.inc'
include 'epos.incsem'
common/photrans/phoele(4),ebeam,noevt
common/photrans/phoele(4),ebeam
double precision pgampr,rgampr
common/cgampr/pgampr(5),rgampr(4)
......@@ -56,13 +56,13 @@ c PHOJET common
& ,ELEM,ELEM2,Q2MN,Q2MX,XIMAX,XIMIN,XIDEL,DELLY
& ,FLUXT,FLUXL,Q2LOW,Y,FFT,FFL,AY,AY2,YY,WGMAX
& ,ECMIN2,ECMAX2,Q22MIN,Q22AVE,Q22AV2,Q22MAX,AN2MIN
& ,AN2MAX,YY2MIN,YY2MAX,EEMIN2
& ,AN2MAX,YY2MIN,YY2MAX,EEMIN2,gamNsig0
& ,Q21MIN,Q21MAX,AN1MIN,AN1MAX,YY1MIN,YY1MAX
SAVE EE1,EE2,PROM2,YMIN2,YMAX2,THMIN2,THMAX2
& ,ELEM,ELEM2,Q2MN,Q2MX,XIMAX,XIMIN,XIDEL,DELLY
& ,FLUXT,FLUXL,Q2LOW,Y,FFT,FFL,AY,AY2,YY,WGMAX
& ,ECMIN2,ECMAX2,Q22MIN,Q22AVE,Q22AV2,Q22MAX,AN2MIN
& ,AN2MAX,YY2MIN,YY2MAX,EEMIN2
& ,AN2MAX,YY2MIN,YY2MAX,EEMIN2,gamNsig0
& ,Q21MIN,Q21MAX,AN1MIN,AN1MAX,YY1MIN,YY1MAX
INTEGER ITRY,ITRW
SAVE ITRY,ITRW
......@@ -70,7 +70,7 @@ c PHOJET common
C LOCAL VARIABLES
DOUBLE PRECISION PINI(5),PFIN(5),GGECM,PFTHE,YQ2,YEFF,Q2LOG,WGH,Q2
& ,WEIGHT,Q2E,E1Y,PHI,COF,SIF,WGY,DAY,P1(5),P2(4)
& ,drangen
& ,drangen,ALLM97,gamNsig
......@@ -134,7 +134,7 @@ C
FLUXT = FLUXT*DELLY
FLUXL = FLUXL*DELLY
WRITE(ifch,'(1X,A,1P2E12.4)')
& 'PHO_GPHERA: integrated flux (trans./long.):',FLUXT,FLUXL
& 'PHOGPHERA: integrated flux (trans./long.):',FLUXT,FLUXL
ENDIF
C
......@@ -159,22 +159,10 @@ C
YY2MAX = 0.D0
ITRY = 0
ITRW = 0
gamNsig0 = 5d0 * ALLM97(Q2LOW,WGMAX)
elseif(imod.eq.1)then !event
if(ish.ge.2.and.noevt.eq.1)then !statistic
C statistics (restore previous value because event was not used (no col)
AY = AY-YY
AY2 = AY2-YY*YY
YY2MIN = YY1MIN
YY2MAX = YY1MAX
Q22MIN = Q21MIN
Q22MAX = Q21MAX
Q22AVE = Q22AVE-Q2
Q22AV2 = Q22AV2-Q2*Q2
AN2MIN = AN1MIN
AN2MAX = AN1MAX
endif
C
C sample y
ITRY = ITRY+1
......@@ -261,6 +249,10 @@ C ECMS cut
& -(P1(2)+P2(2))**2-(P1(3)+P2(3))**2
IF((GGECM.LT.ECMIN2).OR.(GGECM.GT.ECMAX2)) GOTO 175
GGECM = SQRT(GGECM)
C accept A2 and W according to gamma-p cross section (function of F2)
gamNsig=ALLM97(Q2,GGECM)/gamNsig0
if(gamNsig.ge.1d0)print *,'R>1 in DIS',gamNsig
if(drangen(gamNsig).gt.gamNsig)goto 175 !no interaction
C output
engy=sngl(GGECM)