1 c-----------------------------------------------------------------------
3 c-----------------------------------------------------------------------
4 c energy-momentum sharing
5 c-----------------------------------------------------------------------
9 common/cwzero/wzero,wzerox
10 double precision omega,omlog,oma,omb,wab,wba,wmatrix,wzero,nbar
11 *,wzerox,rrr,eps,xprem,xmrem,om1intgck
13 common/col3/ncol,kolpt
16 double precision s,px,py,pomass,plc!,PhiExpo
17 common/emsptl/nppr(npommx,kollmx),npproj(mamx),nptarg(mamx)
18 common/cncl3/iactpr(mamx),iacttg(mamx)
21 dimension ishuff(500),icp(2),ict(2)
22 call utpri('emsaa ',ish,ishini,4)
38 call emsipt !initialize projectile and target
39 call emsigr !initialize grid
43 if((iactpr(1).eq.1.and.iacttg(1).eq.1.and.maproj+matarg.eq.2)
45 &.or.(iactpr(1).eq.1.and.maproj.eq.1.and.maproj+matarg.gt.2)
46 &.or.(iacttg(1).eq.1.and.matarg.eq.1.and.maproj+matarg.gt.2)
47 &.or.(maproj.gt.1.and.matarg.gt.1))then !if not nothing
54 nSprmx=nSprmx+nprmx(k)
59 if(nemsi.le.4.and.iemsi1.eq.1)call xEmsI1(1,0,omlog)
60 if(ish.ge.6)write (ifch,*)'after xEmsI1'
61 if(nemsi.le.4.and.iemsi2.eq.1)call xEmsI2(1,0)
62 if(ish.ge.6)write (ifch,*)'after xEmsI2'
63 if(ish.ge.6)call XPrint('Before Markov:&')
69 if(ish.ge.4)write(ifch,*)'Markov Process'
70 kint=int(max(15.,2.*engy**0.2))
71 if(koll.gt.50)kint=3*kint/int(log(float(koll)))
72 kmcmx=nSprmx*kint !50*kint !100*kint
75 do kmc=1,kmcmx !-----> start Metropolis
80 knprmx=knprmx+nprmx(ik)
81 if(rrr.le.dble(knprmx)/dble(nSprmx))then ! k-th pair
90 n=1+int(rangen()*float(nprmx(k))) ! n-th spot for k-th pair
92 if(idpr(n,k).eq.0)nbar=nbar-1d0
94 xprem=1.d0!xpp(ip)+xppr(n,k) !consistently, it should be 1.
95 xmrem=1.d0!xmt(it)+xmpr(n,k)
97 wzero=wzerox / ( wzerox
98 & +om1intgck(k,xprem,xmrem)*gammaV(k) )
100 if(ish.ge.8)write(ifch,*)'wzero',k,n,wzero,wzerox,gammaV(k)
101 & ,om1intgck(k,xprem,xmrem)
102 if(ish.ge.1.and.100000*(kmc/100000).eq.kmc)
103 & write(ifmt,*)'kmc',kmc,kmcmx
112 if(idpr(n,k).eq.0.and.idx0.eq.0)then
124 if(ish.ge.8)write(ifch,*)'omb',omb,wab,k,n
126 write (ifmt,*)'wab,kmc',wab,omb,kmc,k,n,xpr(n,k),ypr(n,k)
127 & ,xppr(n,k),xmpr(n,k),xpp(ip),xmt(it),ip,it,idpr(n,k)
128 write(ifmt,'(a,i12,d25.15)')'ems,seedf',nrevt+1,seedc
136 if(oma.ge.0.d0.and.oma.le.eps*omb*wba/wab)then
140 omlog=omlog+dlog(omb)
142 elseif(oma.le.1.d-300.or.oma.ne.oma.or.omb.ne.omb)then
143 write (ifmt,*)'oma,kmc',oma,omb,kmc,k,n,xpr(n,k),ypr(n,k)
144 & ,xppr(n,k),xmpr(n,k),idpr(n,k),npr(1,k),xpp(ip),xmt(it),ip,it
145 write(ifmt,'(a,i12,d25.15)')'ems,seedf',nrevt+1,seedc
150 z=sngl(omb/oma*wba/wab)
151 if(ish.ge.8)write(ifch,*)'z,oma',z,oma,wba,k,n
152 if(rangen().gt.z)then
158 omlog=omlog-dlog(oma)+dlog(omb)
168 kplot=int(float(kmc)/float(kmcmx)*100.)
169 if(iemsi1.eq.1)call xEmsI1(1,kplot,omlog)
170 if(iemsi2.eq.1)call xEmsI2(1,kplot)
173 enddo !-----> end Metropolis
176 elseif(iokoll.eq.1)then
195 c --- Plot Pomeron b-distributions ---
197 if(ish.ge.6)call XPrint('After Markov :&')
201 if(iemsb.eq.1)then ! plot
204 if(nprt(k).gt.0)call xEmsB(1,2,k)
208 if(iemsbg.eq.1.and.ievt0.eq.0)then ! plot
215 if(idpr(n,k).ne.0)call xEmsBg(1,idpr(n,k),k)
221 c --- Plot distr of pomeron number ---
224 if(iemspm.eq.1.and.ievt0.eq.0)then
226 call xEmsPm(1,k,nprt(k))
232 c -- Split Enhanced Pomerons and fix their nature ---
236 if(idfpr(n,k).eq.0)call ProPoTy(k,n)
241 c --- Count real interactions ---
250 kolp(ip)=kolp(ip)+nprt(k) !number of cut Pomerons
251 kolt(it)=kolt(it)+nprt(k) !on remnants
252 c kolp(ip)=kolp(ip)+1
253 c kolt(it)=kolt(it)+1
258 if(ish.ge.5)write(ifch,*)'ncol:',ncol
261 c --- Calculate Z (written to zzremn)
272 c --- recalculate Zptn
275 if(irzptn.eq.1)call recalcZPtn
278 c --- fix all variables
281 if(ish.ge.4)write(ifch,*)'fix all variables'
284 if(itpr(k).ge.2)call ProDiSc(k)
288 if(lproj(ip).ne.0)call ProReEx( 1,ip)
291 if(ltarg(it).ne.0)call ProReEx(-1,it)
300 anintdiff=anintdiff+1.
301 if((iep(1).eq.0.and.iet(1).eq.2).or.
302 & (iet(1).eq.0.and.iep(1).eq.2))anintsdif=anintsdif+1.
309 if(aidif.ge.0..and.itpr(k).eq.2)then
310 aidif=aidif+iep(k)+iet(k)
311 elseif(itpr(k).eq.1)then
316 if(aidif.gt.0.)anintdiff=anintdiff+1.
317 if(aidif.eq.2.)anintsdif=anintsdif+1.
318 if(aiine+aidif.gt.0.)anintine=anintine+1.
322 if(ish.ge.6)call XPrint('After fixing:&')
325 c --- Plot MC pomeron number ---
327 if(nemsi.le.4.and.ntry2.eq.0.and.irea.ge.0)then
328 if(iemsi1.eq.1)call xEmsI1(1,100,omlog)
329 if(iemsi2.eq.1)call xEmsI2(1,100)
330 if(iemsi1.eq.1.and.ncol.gt.0)call xEmsI1(2,0,omlog)
331 if(iemsi2.eq.1.and.ncol.gt.0)call xEmsI2(2,0)
332 if((iemsi1.eq.1.or.iemsi2.eq.1).and.ncol.eq.0)nemsi=nemsi-1
336 if(iemsb.eq.1)then ! plot
338 if(itpr(k).eq.0)call xEmsB(1,3,k)
339 if(itpr(k).eq.1)call xEmsB(1,4,k)
340 if(itpr(k).eq.2)call xEmsB(1,5,k)
341 if(itpr(k).ne.0)call xEmsB(1,6,k)
348 if(ncol.eq.0)goto 998
350 c --- Treat Pomerons ---------------------------------------
353 c --- Check minimum mass ---
357 if(xpr(n,k).lt.cumpom**2/engy**2)then
361 nbkpr(nnv,k)=0 !if bckp Pomeron
364 call VirPom(k,nnb,1) !if hard backup exist
365 nbkpr(n,k)=0 !remove it
372 c --- Set String End Type and Pt
379 if(idpr(n,k).gt.0)then
389 if(ntry.ge.200)vpom=.true.
391 if(ish.ge.4)write(ifch,*)'Try again setting string ends for k,n'
405 c --- Check Pomeron mass
409 if(idpr(n,k).ne.0.and.ivpr(n,k).ne.0)then
410 px=xxp1pr(n,k)+xxp2pr(n,k)+xxm1pr(n,k)+xxm2pr(n,k)
411 py=xyp1pr(n,k)+xyp2pr(n,k)+xym1pr(n,k)+xym2pr(n,k)
412 pomass=xpr(n,k)*plc*plc-px*px-py*py
413 if(pomass.le.0.d0)then
416 idfpom=iabs(idfpr(n,k))
418 call VirPom(k,n,3) !call RmPt(k,n)
419 if(nnv.ne.0)then !bckp Pomeron
422 if(nnb.ne.0)then !big Pomeron with bckp one
437 c --- Define String ends for "backup" Pomerons ---
441 if(nvpr(n,k).ne.0)call ProSeX(k,n,iret)
456 c --- Define String ends for "normal" Pomerons ---
460 if(nvpr(n,k).eq.0)call ProSeX(k,n,iret)
480 if(ncol.eq.0)goto 1000
483 if(itpr(k).eq.1)call emswrpom(k,iproj(k),maproj+itarg(k))
487 c --- Treat hard Pomeron
491 if(idpr(n,k).eq.3)then
493 call psahot(k,n,iret)
495 if(nbkpr(n,k).ne.0)then
497 call ProSeX(k,nn,iret2)
500 istptl(nppr(nn,k))=32
505 idfpr(nn,k)=idfpr(n,k)
508 ansff=ansff+1 !counters
514 elseif(nbkpr(n,k).ne.0)then
517 istptl(nppr(nn,k))=32
523 if(nbkpr(n,k).ne.0)then
525 istptl(nppr(nn,k))=32
533 c --- Treat "normal" soft Pomerons ---
537 if(nvpr(n,k).eq.0)then
539 call ProSeF(k,n,iret)
553 c --- Treat Remnants -----------------------------------------
563 c Here and later "kolp(ip).ne.0" replaced by "iep(ip).ne.-1" to count
564 c projectile and target nucleons which are counted in paires but are not used
565 c in collision (no diffractive or inelastic interaction) as slow particles
566 c at the end. Then we can use them in ProRem to give mass to all other nucleons
567 c and avoid energy conservation violation that utrescl can not treat
568 c (and it gives a reasonnable number of grey particles even if distributions
569 c are not really reproduced).
570 c if(kolp(ip).ne.0)call ProCop(ip,ip)
571 if(iep(ip).ne.-1)call ProCop(ip,ip)
574 if(iet(it).ne.-1)call ProCot(it,maproj+it)
575 c if(kolt(it).ne.0)call ProCot(it,maproj+it)
579 c ---- Remnant Masses (ProReM)
582 if(ish.ge.6)call XPrint('Before ProReM:&')
584 call StoRe(1) !Store Remnant configuration
586 nshuffi=maproj+matarg
588 &call utstop('ems: increase dimension for ishuff&')
590 ishuff(ip)=ip !positive for projectile
593 ishuff(maproj+it)=-it !negative for target
598 do while(nshuff.gt.0)
600 if(nshuff.gt.1.and.koll.eq.1.and.iep(1).ne.iet(1))then
601 c to set the mass of diffractive or not excited remnant first
602 c (to avoid unlimited mass of inelastic remants)
603 if(iep(1).ne.1.and.iep(1).lt.iet(1))then
609 indx=1+int(rangen()*float(nshuff))
611 if(ishuff(indx).gt.0)then
613 c if(kolp(ip).ne.0)call ProReM( 1,ip,iret)
614 if(iep(ip).ne.-1)call ProReM( 1,ip,iret)
617 c if(kolt(it).ne.0)call ProReM(-1,it,iret)
618 if(iet(it).ne.-1)call ProReM(-1,it,iret)
622 !----------------------------------------
623 !If there is a problem, try again shuffle (30 times),
624 !if it doesn't work, for pp, try 10 times with the same type
625 !of event and if doesn't work redo event;
626 !for pA redo event ; and for AB (with A or B >10)
627 !continue with some ghosts ...
628 !----------------------------------------
630 if(ish.ge.3)write(ifch,*)'shuffle, try again',ntry
631 call StoRe(-1) !Restore Remnant configuration
634 elseif(maproj+matarg.eq.2.and.ntry2.lt.10)then
636 if(ish.ge.2)write(ifch,*)'ProRem, try again ! ntry=',ntry2
637 if(ish.ge.2)write(ifmt,*)'ProRem, try again ! ntry=',ntry2
639 elseif(maproj.le.10.or.matarg.le.10)then
640 if(ish.ge.2)write(ifch,*)'ProRem, redo event ! ntry=',ntry
641 if(ish.ge.2)write(ifmt,*)'ProRem, redo event ! ntry=',ntry
648 ishuff(indx)=ishuff(nshuff)
654 if(ish.ge.6)call XPrint('After ProReM:&')
657 c --- Write Z into zpaptl for connected strings
662 if(kolp(ip).ne.0)call WriteZZ(1,ip)
665 if(kolt(it).ne.0)call WriteZZ(-1,it)
674 c if(kolp(ip).ne.0)call emswrp(ip,ip)
675 if(iep(ip).ne.-1)call emswrp(ip,ip)
678 c if(kolt(it).ne.0)call emswrt(it,maproj+it)
679 if(iet(it).ne.-1)call emswrt(it,maproj+it)
683 c --- Remnant Flavors (ProReF)
700 if(nprt(ko).gt.0)then
702 if(idpr(np,ko).gt.0)then
703 call xEmsPx(1,sngl(xpr(np,ko)),sngl(ypr(np,ko)),nprt(ko))
714 if(idpr(n,k).eq.3)then
715 je1=min(1,nemispr(1,n,k))
716 je2=min(1,nemispr(2,n,k))
718 call xEmsP2(1,1+idhpr(n,k),jex
721 * ,sngl(xpprbor(n,k)),sngl(xmprbor(n,k))
722 * ,ptprboo(1,n,k),ptprboo(2,n,k) )
732 if(nprt(ko).gt.0)then
734 if(idpr(np,ko).gt.0)then
735 ptp1=sngl(xxp1pr(np,ko)**2+xyp1pr(np,ko)**2)
736 ptp2=sngl(xxp2pr(np,ko)**2+xyp2pr(np,ko)**2)
737 ptm1=sngl(xxm1pr(np,ko)**2+xym1pr(np,ko)**2)
738 ptm2=sngl(xxm2pr(np,ko)**2+xym2pr(np,ko)**2)
739 call xEmsSe(1,sngl(xp1pr(np,ko)),ptp1,1,1)
740 call xEmsSe(1,sngl(xp2pr(np,ko)),ptp2,1,1)
741 call xEmsSe(1,sngl(xm1pr(np,ko)),ptm1,-1,1)
742 call xEmsSe(1,sngl(xm2pr(np,ko)),ptm2,-1,1)
743 call xEmsSe(1,sngl(xp1pr(np,ko)),sngl(xm1pr(np,ko)),1,2)
744 call xEmsSe(1,sngl(xm2pr(np,ko)),sngl(xp2pr(np,ko)),1,2)
752 do i=maproj+matarg+1,nptl
753 if(istptl(iorptl(i)).eq.41)then
754 xpdr=(pptl(4,i)+pptl(3,i))/sngl(plc)
755 xmdr=(pptl(4,i)-pptl(3,i))/sngl(plc)
756 if(ityptl(i).eq.41)call xEmsDr(1,xpdr,xmdr,1)
757 if(ityptl(i).eq.51)call xEmsDr(1,xpdr,xmdr,2)
758 if(ityptl(i).eq.42)call xEmsDr(1,xpdr,xmdr,3)
759 if(ityptl(i).eq.52)call xEmsDr(1,xpdr,xmdr,4)
766 if(kolp(i).gt.0)call xEmsRx(1,1,sngl(xpp(i)),sngl(xmp(i)))
769 if(kolt(j).gt.0)call xEmsRx(1,2,sngl(xmt(j)),sngl(xpt(j)))
778 c write(*,*)'emsaa-iret',iret
779 if(ish.ge.1.and.iret.ne.0)write(ifch,*)'iret not 0 (ems)=> redo'
780 call utprix('emsaa ',ish,ishini,4)
785 c----------------------------------------------------------------------
786 subroutine StoCon(mode,k,n)
787 c----------------------------------------------------------------------
788 c store or restore configuration
789 c mode = 1 (store) or -1 (restore)
790 c k = collision index
792 c----------------------------------------------------------------------
795 include 'epos.incems'
820 elseif(mode.eq.2)then
840 elseif(mode.eq.-1)then
860 elseif(mode.eq.-2)then
881 call utstop('mode should integer from -2 to 2 (without 0)&')
886 c-------------------------------------------------------------------------
887 subroutine RemPom(k,n)
888 c-------------------------------------------------------------------------
890 c-------------------------------------------------------------------------
891 include 'epos.incems'
896 npr(idpr(n,k),k)=npr(idpr(n,k),k)-1 !nr of pomerons
897 nprt(k)=npr(1,k)+npr(3,k)
898 if(idpr(n,k).gt.0)then
899 npp(ip)=npp(ip)-1 !nr of pomerons per proj
900 npt(it)=npt(it)-1 !nr of pomerons per targ
902 xpp(ip)=xpp(ip)+xppr(n,k)
903 xmt(it)=xmt(it)+xmpr(n,k)
915 c-------------------------------------------------------------------------
916 subroutine ProPo(k,n)
917 c-------------------------------------------------------------------------
918 c propose pomeron type = idpr(n,k
919 c-------------------------------------------------------------------------
921 include 'epos.incems'
922 double precision wzero,wzerox
923 common/cwzero/wzero,wzerox
930 if(dble(rangen()).gt.wzero)then
934 c nbr of pomerons per proj
936 c nbr of pomerons per targ
941 npr(idpr(n,k),k)=npr(idpr(n,k),k)+1 !nr of pomerons
942 nprt(k)=npr(1,k)+npr(3,k)
948 c-------------------------------------------------------------------------
949 subroutine ProXY(k,n)
950 c-------------------------------------------------------------------------
951 c propose pomeron x,y
952 c-------------------------------------------------------------------------
955 include 'epos.incems'
956 include 'epos.incsem'
957 double precision xp,xm,om1xprk,om1xmrk,anip,anit,eps
959 parameter (eps=1.d-30)
970 if(idpr(n,k).ne.0)then
973 c because of fom, it's not symetric any more if we choose always xp first
974 c and then xm ... so choose it randomly.
975 if(rangen().lt.0.5)then
976 xp=om1xprk(k,xprem,xmrem,1)
977 xm=om1xmrk(k,xp,xprem,xmrem,1)
979 xm=om1xprk(k,xmrem,xprem,-1)
980 xp=om1xmrk(k,xm,xmrem,xprem,-1)
984 if(xm.gt.eps.and.xp.gt.eps)then
985 ypr(n,k)=0.5D0*dlog(xp/xm)
989 if(ish.ge.1)write(ifmt,*)'Warning in ProXY ',xp,xm
990 npr(idpr(n,k),k)=npr(idpr(n,k),k)-1
992 npr(idpr(n,k),k)=npr(idpr(n,k),k)+1
997 nprt(k)=npr(1,k)+npr(3,k)
998 npp(ip)=npp(ip)-1 !nr of pomerons per proj
999 npt(it)=npt(it)-1 !nr of pomerons per targ
1002 c Update xp and xm of remnant, and change the limit to have big enought mass.
1006 xpp(ip)=xpp(ip)-xppr(n,k)
1007 xppmn(ip)=min(1.d0,anip*xpmn(ip)/xmpmx(ip))
1008 xmt(it)=xmt(it)-xmpr(n,k)
1009 xmtmn(it)=min(1.d0,anit*xtmn(it)/xptmx(it))
1015 c-------------------------------------------------------------------------
1016 double precision function wmatrix(k,n)
1017 c-------------------------------------------------------------------------
1018 c proposal matrix w(a->b), considering pomeron type, x, y
1019 c-------------------------------------------------------------------------
1021 include 'epos.incems'
1022 include 'epos.incsem'
1023 double precision wzero,wzerox,Womegak,xprem,xmrem,om1intgck
1024 common/cwzero/wzero,wzerox
1030 if(idpr(n,k).eq.0)then
1033 xprem=1.d0!xpp(ip)+xppr(n,k)
1034 xmrem=1.d0!xmt(it)+xmpr(n,k)
1035 wmatrix=(1d0-wzero)/om1intgck(k,xprem,xmrem)
1036 * *Womegak(xppr(n,k),xmpr(n,k),xprem,xmrem,k)
1042 c-------------------------------------------------------------------------
1043 double precision function omega(n,k)
1044 c-------------------------------------------------------------------------
1045 c calculates partial omega for spot (k,n)
1046 c-------------------------------------------------------------------------
1049 include 'epos.incems'
1050 include 'epos.incsem'
1051 common/cwzero/wzero,wzerox
1052 double precision wzero,wzerox,eps
1053 parameter(eps=1.d-15)
1054 double precision PhiExpoK,omGamk,xp,xm,fom
1055 double precision plc,s
1057 common/nucl3/phi,bimp
1064 if(xpp(ip).lt.xppmn(ip)+eps.or.xpp(ip).gt.1.d0+eps)goto 1001
1065 if(xmt(it).lt.xmtmn(it)+eps.or.xmt(it).gt.1.d0+eps)goto 1001
1067 omega=xpp(ip)**dble(alplea(iclpro))
1068 & *xmt(it)**dble(alplea(icltar))
1075 c if(idpr(n,k).gt.0)then
1076 c if(nprt(k).le.0)stop'omega: nprt(k) should be positive !!!! '
1081 c bglaub2=sigine/10./pi !10= fm^2 -> mb
1082 c bglaub=sqrt(bglaub2)
1083 c b2x=epscrp*epscrp*bglaub2
1085 c ztgx=epscrw*exp(-b2/2./b2x)*fscra(engy/egyscr)
1086 c zpjx=epscrw*exp(-b2/2./b2x)*fscra(engy/egyscr)
1093 c if(b2.le.bglaub2)nctg=nctg+1
1094 c ztg=ztg+epscrw*exp(-b2/2./b2x)*fscro(engy/egyscr)
1095 c nlpop=nlpop+nprt(kk)
1102 c if(b2.le.bglaub2)ncpj=ncpj+1
1103 c zpj=zpj+epscrw*exp(-b2/2./b2x)*fscro(engy/egyscr)
1104 c nlpot=nlpot+nprt(kk)
1108 ! zpjx+zpj is equal to zparpro(k)
1109 ! ztgx+ztg is equal to zpartar(k)
1110 zprj=zparpro(k) !zsame+zpj
1111 ztgt=zpartar(k) !zsame+ztg
1112 if(idpr(n,k).eq.0)then
1117 c !-------------------------------------------------------------------------
1118 c ! fom : part of Phi regularization; Phi -> Phi^(n) (n = number of Poms)
1119 c ! Phi^(0) relevant for Xsect unchanged, apart of (maybe) normalization (Z)
1120 c !-------------------------------------------------------------------------
1121 omega=omega*omGamk(k,xp,xm)*gammaV(k)*fom(zprj,xm,bk(k))
1122 & *fom(ztgt,xp,bk(k))
1125 omega=omega*PhiExpoK(k,xpp(ip),xmt(it))
1128 if(omega.le.0.d0)goto 1001
1133 if(itarg(kk).ne.it)then
1136 omega=omega*PhiExpoK(kk,xpp(ipl),xmt(itl))
1137 if(omega.le.0.d0)goto 1001
1142 if(iproj(kk).ne.ip)then
1145 omega=omega*PhiExpoK(kk,xpp(ipl),xmt(itl))
1146 if(omega.le.0.d0)goto 1001
1151 if(omega.lt.1.d-100)then
1152 if(ish.ge.6)write(*,*)'omega-exit',omega
1154 elseif(omega.gt.1.d100)then
1155 if(ish.ge.6)write(*,*)'omega-exit',omega
1168 c-------------------------------------------------------------------------
1169 double precision function fom(z,x,b)
1170 c-------------------------------------------------------------------------
1172 double precision x,u,w,z0
1173 !----------------------------------------------------------------
1174 ! part of Phi regularization; Phi -> Phi^(n) (n = number of Poms)
1175 ! Phi^(0) relevant for Xsect unchanged
1176 !----------------------------------------------------------------
1178 if(z.gt.0..and.alpfomi.gt.0.)then
1181 c u=z0*dble(z/z0)**2.
1182 w=u/z0*exp(-dble(b*b/delD(1,iclpro,icltar)))
1184 !---------------------------------------------------
1185 !e=exp(-0.05*u) !analytic function with e(0)=1
1186 !fom=((1-u)+(u+w)*sqrt(x**2+((u-1+e)/(u+w))**2))
1187 ! fom(z=0)=1 fom(x=0)=e fom(x=1)~w
1188 !---------------------------------------------------
1189 fom=1.d0+w*x**betfom
1190 !---------------------------------------------------
1194 c-------------------------------------------------------------------------
1195 subroutine ProPoTy(k,n)
1196 c-------------------------------------------------------------------------
1197 c propose pomeron type
1198 c-------------------------------------------------------------------------
1201 include 'epos.incems'
1202 include 'epos.incsem'
1203 double precision ww,w0,w1,w2,w3,w4,w5,w(0:7),aks,eps
1205 parameter(eps=1.d-10)
1207 dimension nnn(3),kkk(3)
1209 if(idpr(n,k).eq.0)return
1210 if(ish.ge.4)write(ifch,*)'ProPoTy:k,n,idpr,x',k,n,idpr(n,k)
1212 if(idpr(n,k).ne.1)call utstop('ProPoTy: should not happen&')
1230 if(iep(ip).ne.-1)iep(ip)=-1
1231 if(iet(it).ne.-1)iet(it)=-1
1243 call WomTy(w,xh,yp,k)
1246 if(w(0).gt.0.d0)w0=w(0)
1247 if(w(1).gt.0.d0)w1=w(1)
1248 if(w(2).gt.0.d0)w2=w(2)
1249 if(w(3).gt.0.d0)w3=w(3)
1250 if(w(4).gt.0.d0)w4=w(4)
1251 if(w(5).gt.0.d0)w5=w(5)
1253 ww=w0+w1+w2+w3+w4+w5
1254 if(ish.ge.4)write(ifch,*)'ProPoTy:ww,ww_i'
1255 * ,ww,w0/ww*100.d0,w1/ww*100.d0,w2/ww*100.d0
1256 * ,w3/ww*100.d0,w4/ww*100.d0,w5/ww*100.d0
1259 aks=dble(rangen())*ww
1261 if(ww.lt.eps.or.aks.le.w0)then !soft pomeron
1263 if(ish.ge.8)write(ifch,*)'ProPoTy:idpr',idpr(n,k)
1265 elseif(aks.ge.ww-w5)then !diffractive interaction
1267 if(ish.ge.8)write(ifch,*)'ProPoTy:itpr',itpr(k)
1269 npr(0,k)=npr(0,k)+1 !nr of empty cells
1274 if(ish.ge.8)write(ifch,*)'ProPoTy:idpr',idpr(n,k)
1280 if(aks.le.w1)then !gg-pomeron
1282 elseif(aks.le.w1+w2)then !qg-pomeron
1284 elseif(aks.le.w1+w2+w3)then !gq-pomeron
1286 elseif(aks.le.w1+w2+w3+w4)then !qq-pomeron
1289 call utstop('ems-unknown pomeron&')
1294 if(xpr(n,k).gt.0.d0)then
1297 if(idpr(n,k).eq.1)then
1301 if(idpr(n,k).eq.3)then
1309 if(nnn(i).ne.0.and.kkk(i).ne.0.and.cont)then
1311 if(idpr(nnn(i),kkk(i)).eq.3)then
1313 !Backup soft Pomeron if sh not possible later
1320 if(idpr(nn,kb).eq.0)then !empty spot
1325 xpr(nn,kb)=xpr(nb,kb)
1326 ypr(nn,kb)=ypr(nb,kb)
1327 xppr(nn,kb)=xppr(nb,kb)
1328 xmpr(nn,kb)=xmpr(nb,kb)
1329 idfpr(nn,kb)=-idfpr(nb,kb)
1330 bhpr(nn,kb)=bhpr(nb,kb)
1350 if(ish.ge.2)write(ifmt,*)'no empty lattice site, backup lost'
1360 c-------------------------------------------------------------------------
1361 subroutine ProDiSc(k)
1362 c-------------------------------------------------------------------------
1363 c propose diffractive scattering
1364 c-------------------------------------------------------------------------
1367 include 'epos.incsem'
1368 include 'epos.incems'
1369 common/col3/ncol,kolpt
1374 kolp(ip)=kolp(ip)+itpr(k)/2 !number of cut on remnants
1375 kolt(it)=kolt(it)+itpr(k)/2
1381 c-------------------------------------------------------------------------
1382 subroutine ProReEx(ir,ii)
1383 c-------------------------------------------------------------------------
1384 c propose remnant excitation
1385 c for proj (iep) if ir=1 or target (iet) if ir=-1:
1386 c 0 = no, 1 = inel excitation, 2 = diffr excitation
1387 c-------------------------------------------------------------------------
1390 include 'epos.incsem'
1391 include 'epos.incems'
1392 common/cncl3/iactpr(mamx),iacttg(mamx)
1395 if(ir.eq.1)then !proj
1398 if(iactpr(ip).eq.0)stop'ProReEx: should not happen (1)'
1403 if(itpr(kp).eq.1)mine=1
1404 if(itpr(kp).eq.2)mdif=1
1408 if(mine.eq.1)then !inelastic
1409 c if(r.lt.1.-(1.-rexndi(iclpro))**kolp(ip))iep(ip)=1
1410 if(r.lt.1.-(1.-rexndi(iclpro)))iep(ip)=1
1411 elseif(mdif.eq.1)then !diffr
1412 if(r.lt.1.-(1.-rexdif(iclpro))**kolp(ip))iep(ip)=2
1413 c if(r.lt.1.-(1.-rexdif(iclpro)))iep(ip)=2
1415 elseif(ir.eq.-1)then !targ
1418 if(iacttg(it).eq.0)stop'ProReEx: should not happen (2)'
1423 if(itpr(kt).eq.1)mine=1
1424 if(itpr(kt).eq.2)mdif=1
1428 if(mine.eq.1)then !inelastic
1429 c if(r.lt.1.-(1.-rexndi(icltar))**kolt(it))iet(it)=1
1430 if(r.lt.1.-(1.-rexndi(icltar)))iet(it)=1
1431 elseif(mdif.eq.1)then !diffr
1432 if(r.lt.1.-(1.-rexdif(icltar))**kolt(it))iet(it)=2
1433 c if(r.lt.1.-(1.-rexdif(icltar)))iet(it)=2
1441 c-------------------------------------------------------------------------
1442 subroutine ProDiPt(k)
1443 c-------------------------------------------------------------------------
1444 c propose transverse momentum for diffractive interaction
1445 c-------------------------------------------------------------------------
1447 include 'epos.incems'
1449 double precision xxe,xye,p5sqpr,p5sqtg
1450 double precision plc,s
1457 c generate p_t for diffractive
1459 if(ptdiff.ne.0.)then
1460 if(itpr(k).eq.2)then
1464 pt=ranpt()*ptd !sqrt(-alog(r)/ad)
1465 elseif(itpr(k).eq.0)then !pt for non-wounded nucleon (usefull in ProRem to avoid problem in utrescl)
1474 xxe=dble(pt*cos(phi))
1475 xye=dble(pt*sin(phi))
1481 c update remnant p_t
1483 10 xxp(ip)=xxp(ip)-xxe
1488 if(iep(ip).eq.6.and.iet(it).eq.6)then !to simulate the fact that originally
1489 if(itpr(k).eq.3)then
1490 call StoCon(-k,k,1) !to fixe mass of corresponding remnants
1491 xpp(ip)=xpp(ip)-xppr(1,k)
1492 xpt(it)=xpt(it)+xppr(1,k)
1493 xmt(it)=xmt(it)-xmpr(1,k)
1494 xmp(ip)=xmp(ip)+xmpr(1,k)
1501 p5sqpr=xpp(ip)*plc*xmp(ip)*plc-dble(amproj*amproj)
1502 p5sqtg=xpt(it)*plc*xmt(it)*plc-dble(amtarg*amtarg)
1506 pt=ranptcut(ptsemx)*ptsend**2
1507 if(ntry.lt.100.and.(p5sqpr-dble(pt*pt).lt.0.d0
1508 & .or.p5sqtg-dble(pt*pt).lt.0.d0))then
1511 pt=ranptcut(ptsemx)*ptsendi
1513 xxe=dble(pt*cos(phi))
1514 xye=dble(pt*sin(phi))
1524 c-------------------------------------------------------------------------
1525 subroutine ProSePt(k,n)
1526 c-------------------------------------------------------------------------
1527 c propose transverse momentum for string ends
1528 c-------------------------------------------------------------------------
1531 include 'epos.incems'
1533 if(ivpr(n,k).eq.2)return !Backup Pomeron
1537 amk0=(qmass(1)+qmass(2)+qmass(3)) !mass for mt distribution but spoil <pt> of anti-proton at low energy
1539 ctp060829 id=idpr(n,k)
1541 ctp060829 if(id.eq.3)ih=1
1543 c generate p_t for string ends (proj)
1553 c if(isplit.eq.1)then
1554 c if(lproj(ip).ge.1)then
1557 c if(itpr(kpair).eq.1)then
1558 c zz=zz+zparpro(kpair)
1565 c if(iep(ip).eq.0)ptsef=ptsef*ptsendi
1567 if(iep(ip).eq.0)ptsendx=ptsendi
1570 if(idp1pr(n,k).gt.0)then
1571 if(idp1pr(n,k).eq.4.or.idp1pr(n,k).eq.5)then !diquarks
1572 c pt=ranptd()*ptsendx
1573 c if(iep(ip).eq.0)then
1574 c pt=ranpt()*ptsendy
1576 c pt=ranptd()*ptsendy
1578 pt=ranptcut(ptsef)*ptsendy
1579 amk1=amk0!+qmass(0)*0.5 !mass for mt distribution with bounding energy for diquark
1581 c pt=ranptd()*ptsendx
1582 c if(iep(ip).eq.0)then
1583 c pt=ranpt()*ptsendx
1585 c pt=ranptd()*ptsendx
1587 pt=ranptcut(ptsef)*ptsendx
1590 pt=sqrt(pt*pt+2.*pt*amk1) !sample mt-m0 instead of pt ...
1592 xxp1pr(n,k)=dble(pt*cos(phi))
1593 xyp1pr(n,k)=dble(pt*sin(phi))
1598 if(idp2pr(n,k).gt.0)then
1599 if(idp2pr(n,k).eq.4.or.idp2pr(n,k).eq.5)then
1600 c pt=ranptd()*ptsendy
1601 c if(iep(ip).eq.0)then
1602 c pt=ranpt()*ptsendy
1604 c pt=ranptd()*ptsendy
1606 pt=ranptcut(ptsef)*ptsendy
1607 amk1=amk0!+qmass(0)*0.5 !mass for mt distribution with bounding energy for diquark
1609 c pt=ranptd()*ptsendx
1610 c if(iep(ip).eq.0)then
1611 c pt=ranpt()*ptsendx
1613 c pt=ranptd()*ptsendx
1615 pt=ranptcut(ptsef)*ptsendx
1618 pt=sqrt(pt*pt+2.*pt*amk1) !sample mt-m0 instead of pt ...
1620 xxp2pr(n,k)=dble(pt*cos(phi))
1621 xyp2pr(n,k)=dble(pt*sin(phi))
1626 c generate p_t for string ends (targ)
1637 c if(isplit.eq.1)then
1638 c if(ltarg(it).ge.1)then
1641 c if(itpr(kpair).eq.1)then
1642 c zz=zz+zpartar(kpair)
1649 c if(iet(it).eq.0)ptsef=ptsef*ptsendi
1651 if(iet(it).eq.0)ptsendx=ptsendi
1654 if(idm1pr(n,k).gt.0)then
1655 if(idm1pr(n,k).eq.4.or.idm1pr(n,k).eq.5)then
1656 c pt=ranptd()*ptsendy
1657 c if(iet(it).eq.0)then
1658 c pt=ranpt()*ptsendy
1660 c pt=ranptd()*ptsendy
1662 pt=ranptcut(ptsef)*ptsendy
1663 amk1=amk0!+qmass(0)*0.5 !mass for mt distribution with bounding energy for diquark
1665 c pt=ranptd()*ptsendx
1666 c if(iet(it).eq.0)then
1667 c pt=ranpt()*ptsendx
1669 c pt=ranptd()*ptsendx
1671 pt=ranptcut(ptsef)*ptsendx
1674 pt=sqrt(pt*pt+2.*pt*amk1) !sample mt-m0 instead of pt ...
1676 xxm1pr(n,k)=dble(pt*cos(phi))
1677 xym1pr(n,k)=dble(pt*sin(phi))
1682 if(idm2pr(n,k).gt.0)then
1683 if(idm2pr(n,k).eq.4.or.idm2pr(n,k).eq.5)then
1684 c pt=ranptd()*ptsendy
1685 c if(iet(it).eq.0)then
1686 c pt=ranpt()*ptsendy
1688 c pt=ranptd()*ptsendy
1690 pt=ranptcut(ptsef)*ptsendy
1691 amk1=amk0!+qmass(0)*0.5 !mass for mt distribution with bounding energy for diquark
1693 c pt=ranptd()*ptsendx
1694 c if(iet(it).eq.0)then
1695 c pt=ranpt()*ptsendx
1697 c pt=ranptd()*ptsendx
1699 pt=ranptcut(ptsef)*ptsendx
1702 pt=sqrt(pt*pt+2.*pt*amk1) !sample mt-m0 instead of pt ...
1704 xxm2pr(n,k)=dble(pt*cos(phi))
1705 xym2pr(n,k)=dble(pt*sin(phi))
1711 c update backup soft pomeron p_t if exist
1713 if(nbkpr(n,k).ne.0)then
1715 xxp1pr(nn,k)=xxp1pr(n,k)
1716 xyp1pr(nn,k)=xyp1pr(n,k)
1717 xxp2pr(nn,k)=xxp2pr(n,k)
1718 xyp2pr(nn,k)=xyp2pr(n,k)
1719 xxm1pr(nn,k)=xxm1pr(n,k)
1720 xym1pr(nn,k)=xym1pr(n,k)
1721 xxm2pr(nn,k)=xxm2pr(n,k)
1722 xym2pr(nn,k)=xym2pr(n,k)
1727 c update remnant p_t (pomeron)
1728 xxp(ip)=xxp(ip)-xxp1pr(n,k)-xxp2pr(n,k)
1729 xyp(ip)=xyp(ip)-xyp1pr(n,k)-xyp2pr(n,k)
1730 xxt(it)=xxt(it)-xxm1pr(n,k)-xxm2pr(n,k)
1731 xyt(it)=xyt(it)-xym1pr(n,k)-xym2pr(n,k)
1734 write(ifch,*) 'ProSePt'
1735 write(ifch,'(4i14/4d14.3/4d14.3/)')
1736 * idp1pr(n,k),idp2pr(n,k),idm1pr(n,k),idm2pr(n,k)
1737 *,xxp1pr(n,k),xxp2pr(n,k),xxm1pr(n,k),xxm2pr(n,k)
1738 *,xyp1pr(n,k),xyp2pr(n,k),xym1pr(n,k),xym2pr(n,k)
1743 c-----------------------------------------------------------------------
1744 subroutine ProSeX(k,n,iret)
1745 c-----------------------------------------------------------------------
1746 c calculates x of string ends
1747 c-----------------------------------------------------------------------
1750 include 'epos.incems'
1752 double precision s,plc
1753 common/cems10/a(0:ntypmx),b(0:ntypmx),d(0:ntypmx)
1754 double precision a,b,d
1755 *,xp,xm,ap1,ap2,am1,am2,aamin1,aamin2,u
1760 if(itpr(k).ne.1)return
1761 if(idpr(n,k).ne.1.or.ivpr(n,k).eq.0)return
1763 if(idp1pr(n,k).eq.0.and.idp2pr(n,k).eq.0
1764 * .and.idm1pr(n,k).eq.0.and.idm2pr(n,k).eq.0)
1765 *call utstop('no Pomeron in ProSex&')
1773 aamin1=ammn(idp1pr(n,k)+idm2pr(n,k))
1774 aamin2=ammn(idp2pr(n,k)+idm1pr(n,k))
1775 xmn1=(aamin1**2+(xxp1pr(n,k)+xxm2pr(n,k))**2
1776 & +(xyp1pr(n,k)+xym2pr(n,k))**2)/s
1777 xmn2=(aamin2**2+(xxp2pr(n,k)+xxm1pr(n,k))**2
1778 & +(xyp2pr(n,k)+xym1pr(n,k))**2)/s
1784 if(ish.ge.5)write(ifch,*)'Problem in ProSex(k,n)',k,n
1788 1 u=dble(rangen())**(1d0/(1d0+ap1))
1789 if(dble(rangen()).gt.(1d0-u)**ap2)goto1
1792 2 u=dble(rangen())**(1d0/(1d0+am1))
1793 if(dble(rangen()).gt.(1d0-u)**am2)goto2
1797 if(xp1pr(n,k)*xm2pr(n,k).lt.xmn1)then
1799 c fc=xp1pr(n,k)*xm2pr(n,k)/xmn1 !avoid virpom
1800 c if(fc.eq.0.)goto 999
1801 c xp1pr(n,k)=xp1pr(n,k)/sqrt(fc)
1802 c xm2pr(n,k)=xm2pr(n,k)/sqrt(fc)
1804 if(xp2pr(n,k)*xm1pr(n,k).lt.xmn2)then
1806 c fc=xp2pr(n,k)*xm1pr(n,k)/xmn2 !avoid virpom
1807 c if(fc.eq.0.)goto 999
1808 c xp2pr(n,k)=xp2pr(n,k)/sqrt(fc)
1809 c xm1pr(n,k)=xm1pr(n,k)/sqrt(fc)
1813 write(ifch,*) 'ProSeX'
1814 write(ifch,'(2d28.3,i8)') xp,xm,ntry
1815 write(ifch,'(4d14.3)')xp1pr(n,k),xp2pr(n,k),xm1pr(n,k),xm2pr(n,k)
1816 write(ifch,'(4d14.3/)')xp1pr(n,k)*xm2pr(n,k)
1817 * ,xp2pr(n,k)*xm1pr(n,k), xmn1, xmn2
1821 c-------------------------------------------------------------------------
1822 subroutine RmPt(k,n)
1823 c-------------------------------------------------------------------------
1824 c remove pt from pomeron
1825 c-------------------------------------------------------------------------
1827 include 'epos.incems'
1830 xxp(ip)=xxp(ip)+xxp1pr(n,k)+xxp2pr(n,k)
1831 xyp(ip)=xyp(ip)+xyp1pr(n,k)+xyp2pr(n,k)
1832 xxt(it)=xxt(it)+xxm1pr(n,k)+xxm2pr(n,k)
1833 xyt(it)=xyt(it)+xym1pr(n,k)+xym2pr(n,k)
1852 c-------------------------------------------------------------------------
1853 subroutine VirPom(k,n,id)
1854 c-------------------------------------------------------------------------
1855 c create virtual pomeron
1856 c virtual pomeron: ivpr(n,k)=0, otherwise ivpr(n,k)=1
1857 c-------------------------------------------------------------------------
1860 include 'epos.incems'
1861 common/col3/ncol,kolpt
1862 common/emsptl/nppr(npommx,kollmx),npproj(mamx),nptarg(mamx)
1863 double precision plc,s
1868 call utpri('VirPom',ish,ishini,3)
1870 if(idpr(n,k).eq.0)return
1879 c print *,' ',id,' ',nvir
1882 write(ifch,*)"virpom ",id," (n,k)",n,k,nnb,nnv,nppr(n,k)
1883 if(ish.ge.5)write(ifch,*)"remnant in",xpp(ip),xmt(it)
1889 if(idpr(nn,kk).eq.0)then
1897 if(idpr(nn,kk).eq.0)then
1903 if(nbkpr(n,k).eq.0.and.nvpr(n,k).eq.0)then !normal Pomeron
1908 npr(idpr(n,k),k)=npr(idpr(n,k),k)-1
1909 nprt(k)=npr(1,k)+npr(3,k)
1911 if(idpr(n,k).eq.1)ansff=ansff-1
1912 if(idpr(n,k).eq.3)anshf=anshf-1
1915 xpp(ip)=xpp(ip)+xppr(n,k)
1916 xmt(it)=xmt(it)+xmpr(n,k)
1917 xxp(ip)=xxp(ip)+xxp1pr(n,k)+xxp2pr(n,k)
1918 xyp(ip)=xyp(ip)+xyp1pr(n,k)+xyp2pr(n,k)
1919 xxt(it)=xxt(it)+xxm1pr(n,k)+xxm2pr(n,k)
1920 xyt(it)=xyt(it)+xym1pr(n,k)+xym2pr(n,k)
1922 if(itpr(k).eq.1.and.nprt(k).eq.0)then !no more Pomeron on this pair
1923 if(kolp(ip).eq.0)then
1924 kolp(ip)=1 !excite nucleon (remnant get pt in ProDiPt)
1925 iep(ip)=6 !with inel mass and inverted string
1927 if(kolt(it).eq.0)then
1931 if(koll.le.2.and.iep(ip).eq.6.and.iet(it).eq.6)then !for small systems we can store lost informations
1932 itpr(k)=3 !to use it to define remnant mass
1933 call StoCon(k,k,n) !store information on lost Pomeron
1966 if(ish.ge.5)write(ifch,*)"remnant out",xpp(ip),xmt(it)
1968 call utprix('VirPom',ish,ishini,3)
1972 c-----------------------------------------------------------------------
1973 subroutine StoRe(imod)
1974 c-----------------------------------------------------------------------
1975 c Store Remnant configuration (imod=1) before shuffle to restore the
1976 c initial configuration (imod=-1) in case of problem.
1977 c-----------------------------------------------------------------------
1980 include 'epos.incems'
1984 c initialize projectile
2000 elseif(imod.eq.-1)then
2002 c restore projectile
2020 call utstop('Do not know what to do in StoRe.&')
2027 c-----------------------------------------------------------------------
2028 subroutine CalcZZ(ir,m)
2029 c-----------------------------------------------------------------------
2030 C Calculates zz for remnant m for proj (ir=1) or target (ir=-1)
2031 c writes it to zzremn(m, 1 or 2)
2032 c-----------------------------------------------------------------------
2034 include 'epos.incems'
2036 if(kolp(m).eq.0)then
2040 elseif(ir.eq.-1)then
2041 if(kolt(m).eq.0)then
2049 if(lproj(m).ge.1)then
2052 if(itpr(kpair).eq.1)then
2053 zz=zz+zparpro(kpair)
2058 elseif(ir.eq.-1)then
2060 if(ltarg(m).ge.1)then
2063 if(itpr(kpair).eq.1)then
2064 zz=zz+zpartar(kpair)
2070 stop'CalcZZ: invalid option. '
2073 if(ir.eq.1) zzremn(m,1)=0
2074 if(ir.eq.-1)zzremn(m,2)=0
2078 c-----------------------------------------------------------------------
2079 subroutine WriteZZ(ir,irem)
2080 c-----------------------------------------------------------------------
2081 c Write Z into zpaptl(K) for connected strings
2082 c K is the index for the string end
2083 c on the corresponding remnant side
2084 c-----------------------------------------------------------------------
2087 include 'epos.incems'
2088 common/emsptl/nppr(npommx,kollmx),npproj(mamx),nptarg(mamx)
2092 elseif(ir.eq.-1)then
2096 do li=1,lremn(irem,jrem)
2097 kkk=kremn(irem,li,jrem)
2099 if(idpr(n,kkk).ne.0)then
2101 c write(ifch,*)'remn',irem,' (',jrem,' ) pom',npom
2102 c & ,' ',zzremn(irem,jrem)
2104 do is=ifrptl(1,npom),ifrptl(2,npom)
2106 if(idptl(is).ne.9)ie=ie+1
2110 if(ir.eq. 1)zpaptl(is1)=zzremn(irem,jrem)
2111 if(ir.eq.-1)zpaptl(is2)=zzremn(irem,jrem)
2113 c write(ifch,*)' ',isi,idptl(isi),zpaptl(isi)
2123 c-----------------------------------------------------------------------
2124 subroutine ProReM(ir,irem,iret)
2125 c-----------------------------------------------------------------------
2126 c propose remnant mass of remnant irem in case of proj (ir=1)
2129 c iret : input : if iret=10 force to give mass even if no more energy,
2130 c when input not 10 : output = error if 1
2131 c-----------------------------------------------------------------------
2134 include 'epos.incems'
2135 double precision rr,xxx,xmin,xmax,msmin,xmmin,xpt2rem,xtest0
2136 double precision at,alp,xi,xii,eps,sx,xmin0,xtest(mamx),fxtest
2137 parameter(eps=1.d-20)
2138 common/cemsr5/at(0:1,0:6)
2139 double precision plc,s,p5sq,aremn,aremnex
2145 call utpri('ProReM',ish,ishini,5)
2155 c uncomment the following two lines to force the excitation
2157 ccc force=.true. !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
2158 ccc ntrymx=1 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
2160 c initial definitions
2173 !idx=isign(iabs(idproj)/10*10+1,idproj)
2174 !call idmass(idx,amremn)
2176 ctp kolzi=kolp(irem)
2178 msmin=dble(amremn*amremn)
2179 if(iep(irem).eq.6)goto 678
2180 elseif(ir.eq.-1)then
2189 !idx=isign(iabs(idtarg)/10*10+1,idtarg)
2190 !call idmass(idx,amremn)
2192 ctp kolzi=kolt(irem)
2194 msmin=dble(amremn*amremn)
2195 if(iet(irem).eq.6)goto 678
2198 ctp noevt replace noxevt
2199 ctp if iez=0, 5% energy violation allowed to give mass to the other side
2200 ctp noxevt=0 !?????? otherwise, energy is strongly not conserved
2202 ctp if(iez(i,jremo).gt.0)noxevt=noxevt+1
2209 if(ntry.gt.ntrymx)then
2211 call utmsg('ProReM')
2212 write(ifch,*)'Remnant mass assignment not possible (ntry)'
2213 if(force)write(ifch,*)'Ignore p4 conservation'
2230 if(xpz(irem,jrem).le.0.d0)then
2231 write(ifch,*)'ProRem ipp',xpz(irem,jrem),irem,lremn(irem,jrem)
2232 do li=1,lremn(irem,jrem)
2233 kkk=kremn(irem,li,jrem)
2234 write(ifch,*)'kkk',kkk
2236 call XPrint('ProRem pro:&')
2237 call utstop('Big problem in ProRem pro!&')
2240 c xtest = xminus-max, corresponding mostly to a remnant mass 0.2
2247 ctp if(xmz(j,jremo).gt.eps.and.iez(j,jrem).gt.0)then !xmz(,jremo)=xplus
2248 ctp060824 if(xmz(j,jremo).gt.eps.and.iez(j,jrem).ge.0)then !xmz(,jremo)=xplus
2249 c if(iez(j,jremo).gt.0.or.koll.eq.1)then !xmz(,jremo)=xplus
2250 if(xmz(j,jremo).gt.eps)then !xmz(,jremo)=xplus
2252 xmmin=xzos(j,jremo)/xmz(j,jremo)
2256 xtest(j)=xpz(j,jremo)-xmmin !maximal momentum available
2257 !this term is very important for non excited remnants in pp, it changes the xf
2258 ! distribution of proton and the multiplicity at low energy. Fxtest should not
2259 ! be to close to 0. otherwise it makes a step in xf distribution of p at
2260 ! 1-fxtest but if fxtest=1, multiplicity at low energy is too high ...
2262 if(xtest(j).gt.0d0)then
2263 xtest(j)=min(xtest(j),fxtest/xpz(irem,jrem))
2265 xtest(j)=min(1.d0,fxtest/xpz(irem,jrem))
2269 c xtest(j)=0.01d0 !maximal momentum available for non exited state
2271 xtest0=max(xtest0,xtest(j))
2272 c print *,iep(1),iet(1),iez(irem,jrem),xtest(j),xpz(j,jremo),xmmin
2273 c & ,xzos(j,jremo),xmz(j,jremo)
2275 ctp060824 if(.not.cont)xtest=min(1.d0,0.2d0/xpz(irem,jrem))
2283 icrmn(1)=icremn(1,irem,jrem)
2284 icrmn(2)=icremn(2,irem,jrem)
2285 xpt2rem=xxz(irem,jrem)**2d0+xyz(irem,jrem)**2d0
2287 c fremnux (+) or fremnux2 (-) ?
2289 aremn=dble(fremnux2(icrmn)) !dble(max(amremn,fremnux2(icrmn)))
2290 c if(iez(j,jrem).eq.2)then
2293 c aremnex=max(amzmn(idz(irem,jrem),jrem) !makes remnant to heavy at low energy
2294 c & +amemn(idz(irem,jrem),iez(irem,jrem))
2295 c & ,dble(fremnux(icrmn)))
2296 aremnex=aremn+amemn(idz(irem,jrem),iez(irem,jrem))
2301 c xmin0=1.05*(aremn**2d0+xxz(irem,jrem)**2d0+xyz(irem,jrem)**2d0)/sx
2302 c xmin=1.1*(aremnex**2d0+xxz(irem,jrem)**2d0+xyz(irem,jrem)**2d0)/sx
2303 xmin0=1.01d0*(aremn**2d0+xpt2rem)/sx
2304 xmin=1.01d0*(aremnex**2d0+xpt2rem)/sx
2305 xmax=min(1.d6/s,xtest0) !to avoid ultra high mass remnants
2306 c for diffractive remnant, mass should never exceed 5% of the proj or targ energy
2307 if(iez(irem,jrem).eq.1)then
2308 xmax=min(xmax,max(dble(xmaxremn),xmin))
2309 elseif(iez(irem,jrem).eq.2)then
2310 xmax=min(xmax,max(dble(xmaxdiff),xmin))
2312 if(koll.eq.1)xmax=min(xmax,xpz(iremo1,jremo))
2313 if(xmin.ge.xmax-eps)then
2315 if(koll.ne.1)xmax=1.d0
2316 if(xmin.ge.xmax-eps)then
2327 if(iez(irem,jrem).gt.0)then
2328 c xmin=xmin-xpt2rem/sx !no pt
2329 c xmax=xmax-xpt2rem/sx !no pt
2330 alp=at(idz(irem,jrem),iez(irem,jrem))
2331 c print *,'alp',iez(irem,jrem),alp,xmin,xmax
2332 if(dabs(alp-1.d0).lt.eps)then
2333 xxx=xmax**rr*xmin**(1d0-rr)
2335 xxx=(rr*xmax**(1d0-alp)+(1d0-rr)*xmin**(1d0-alp))
2338 c xxx=xxx+xpt2rem/sx !no pt
2340 c xmin=dble(amremn)**2d0/sx !no pt
2341 c xxx=xmin+xpt2rem/sx !no pt
2342 xmin=(dble(amremn)**2d0+xpt2rem)/sx
2344 if(xmin.gt.xmax+eps)then
2345 if(ish.ge.6)write(ifch,*)'xmin>xmax for proj not possible (2)'
2353 !to have nice diffractive pic, do not allow too much fluctuation
2356 xmin0=min(0.99d0,1d0-fxtest*dble(1.-rangen()))*xxx
2357 c xmin0=dble(0.9+0.09*rangen())*xxx
2359 xzos(irem,jrem)=xmin0*xpz(irem,jrem)
2361 c msmin=xmin*sx+xpt2rem !no pt
2363 c partition xminus between nucleons of the other side
2366 iimax=noevt !number of opposite side participants
2367 ctp iimax=noxevt !number of opposite side participants
2369 iro=int(rangen()*masso)+1 ! choose ramdomly a nucleon to start
2373 cont=iez(iro,jremo).lt.0.or.xme(iro).lt.-0.99
2376 if(iro.gt.masso)iro=iro-masso
2379 cont=iez(iro,jremo).lt.0.or.xme(iro).lt.-0.99
2383 xi=xii*dble(rangen())**(1.d0/dble(ii-1))
2387 xme(iro)=xxx*(xii-xi)
2389 xmmin=xzos(iro,jremo)
2390 if(xmz(iro,jremo).gt.eps)then
2391 xmmin=xmmin/xmz(iro,jremo)
2392 elseif(koll.eq.1.and.xtest(iro).gt.eps)then
2393 xmmin=xmmin/min(xpz(irem,jrem),xtest(iro))
2394 elseif(xtest(iro).gt.eps)then
2395 xmmin=xmmin/xtest(iro)
2397 if((xpz(iro,jremo)-xme(iro)).lt.xmmin)then
2398 c write(ifch,*)' skip ',cremn,' ',ii,iimax,ntry,xxx
2399 c & ,xpz(iro,jremo)-xme(iro),xmmin
2404 c write(ifch,*)' ok ',cremn,' ',ii,iimax,ntry,xme(iro)/xxx
2407 if(iro.gt.masso)iro=iro-masso
2412 c check xmz(irem,jrem)
2416 678 p5sq=xpz(irem,jrem)*plc*xmz(irem,jrem)*plc
2417 c write(ifch,*)'final mass',p5sq,msmin,xpz(irem,jrem),xmz(irem,jrem)
2419 if(p5sq.lt.msmin)then
2421 call utmsg('ProReM')
2422 write(ifch,*)'Remnant mass assignment not possible (M<Mmin)!'
2423 if(force)write(ifch,*)'Ignore p4 conservation'
2429 elseif(xpz(irem,jrem).gt.0.d0)then
2430 xmz(irem,jrem)=msmin/(plc*plc*xpz(irem,jrem))
2437 if(xme(iro).gt.0.d0)then
2438 xpz(iro,jremo)=xpz(iro,jremo)-xme(iro) !xpz(,jremo)=xminus
2444 call utprix('ProReM',ish,ishini,5)
2448 c-----------------------------------------------------------------------
2449 subroutine ProSeTy(k,n)
2450 c-----------------------------------------------------------------------
2451 c creates proposal for string ends, idp., idm.
2452 c updates quark counters
2453 c-----------------------------------------------------------------------
2455 include 'epos.incems'
2457 common/ems6/ivp0,iap0,idp0,isp0,ivt0,iat0,idt0,ist0
2458 double precision pes,xfqp,xfqt !so01
2459 parameter(eps=1.e-6)
2460 common/ems9/xfqp(0:9),xfqt(0:9)
2461 common/emsx3/pes(0:3,0:6)
2463 if(idpr(n,k).eq.2)stop'no Reggeons any more'
2468 if(idpr(n,k).eq.3)then
2481 if(idhpr(n,k).eq.3)then !so01
2488 elseif(idhpr(n,k).eq.2)then
2494 elseif(idhpr(n,k).eq.1)then
2500 elseif(idhpr(n,k).eq.0)then
2506 call utstop('ProSeTy-idhpr????&')
2510 elseif(idpr(n,k).eq.1)then
2514 if(iabs(idfpr(n,k)).eq.1)then
2518 if(ntry.gt.10)call utstop('something goes wrong in sr ProSeTy&')
2520 pvs=wgtval ! *ivp(ip)
2521 psa=wgtval ! *iap(ip)
2523 psvv=wgtqqq(iclpro) ! *ivp(ip)*(ivp(ip)-1)/2.
2524 paas=wgtqqq(iclpro) ! *iap(ip)*(iap(ip)-1)/2.
2525 su=pss+pvs+psa+pdd+psvv+paas
2533 if(r.gt.(pssp+pvsp+psap+pddp+psvvp).and.paasp.gt.eps)then
2538 elseif(r.gt.(pssp+pvsp+psap+pddp).and.psvvp.gt.eps)then
2543 elseif(r.gt.(pssp+pvsp+psap).and.pddp.gt.eps)then
2547 elseif(r.gt.(pssp+pvsp).and.psap.gt.eps)then
2552 elseif(r.gt.pssp.and.pvsp.gt.eps)then
2557 elseif(pssp.gt.eps)then
2573 if(iabs(idfpr(n,k)).eq.1)then
2578 if(ntry.gt.10)call utstop('something goes wrong in sr ProSeTy&')
2580 pvs=wgtval ! *ivt(it)
2581 psa=wgtval ! *iat(it)
2583 psvv=wgtqqq(icltar) ! *ivt(it)*(ivt(it)-1)/2.
2584 paas=wgtqqq(icltar) ! *iat(it)*(iat(it)-1)/2.
2585 su=pss+pvs+psa+pdd+psvv+paas
2593 if(r.gt.(psst+pvst+psat+pddt+psvvt).and.paast.gt.eps)then
2598 elseif(r.gt.(psst+pvst+psat+pddt).and.psvvt.gt.eps)then
2603 elseif(r.gt.(psst+pvst+psat).and.pddt.gt.eps)then
2607 elseif(r.gt.(psst+pvst).and.psat.gt.eps)then
2612 elseif(r.gt.psst.and.pvst.gt.eps)then
2617 elseif(psst.gt.eps)then
2631 elseif(idpr(n,k).eq.0)then
2641 write(ifch,'(a,2(6(f3.2,1x),2x),$)')'ProSeTy ',
2642 * pssp,pvsp,psap,pddp,psvvp,paasp, psst,pvst,psat,pddt,psvvt,paast
2643 write(ifch,'(2x,3i3,2x,2(i2,1x,2i2,1x,2i2,2x))')idpr(n,k),n,k
2644 * ,idsppr(n,k),idp1pr(n,k),idp2pr(n,k),ivp(ip),iap(ip)
2645 * ,idstpr(n,k),idm1pr(n,k),idm2pr(n,k),ivt(it),iat(it)
2651 c-----------------------------------------------------------------------
2652 subroutine ProSeF(k,n,iret)
2653 c-----------------------------------------------------------------------
2654 c starting from string properties as already determined in EMS,
2655 c one determines string end flavors
2656 c by checking compatibility with remnant masses.
2657 c strings are written to /cems/ and then to /cptl/
2658 c remnant ic is updated (icproj,ictarg)
2659 c------------------------------------------------------------------------
2662 include 'epos.incems'
2664 double precision plc,s,pstg,pend
2666 common/cems/pstg(5,2),pend(4,4),idend(4)
2667 common/emsptl/nppr(npommx,kollmx),npproj(mamx),nptarg(mamx)
2668 integer icp(2),ict(2),ic(2),icp1(2),icp2(2),icm1(2),icm2(2)
2670 integer jcp1(nflav,2),jcp2(nflav,2),jcm1(nflav,2),jcm2(nflav,2)
2671 common/col3/ncol,kolpt /cfacmss/facmss /cts/its
2673 call utpri('ProSeF',ish,ishini,6)
2681 if(itpr(k).ne.1)return
2686 if(idpr(n,k).eq.0.or.ivpr(n,k).eq.0)return
2687 if(idpr(n,k).eq.2)stop'Reggeon'
2688 if(idpr(n,k).eq.3)goto 1000
2690 write(ifch,*)'soft Pomeron'
2691 write(ifch,*)'k:',k,' n:',n,' ip:',ip,' it:',it
2697 pend(1,1)=xxp1pr(n,k)
2698 pend(2,1)=xyp1pr(n,k)
2699 pend(3,1)=xp1pr(n,k)*plc/2d0
2700 pend(4,1)=dsqrt(pend(1,1)**2+pend(2,1)**2+pend(3,1)**2)
2701 pend(1,2)=xxp2pr(n,k)
2702 pend(2,2)=xyp2pr(n,k)
2703 pend(3,2)=xp2pr(n,k)*plc/2d0
2704 pend(4,2)=dsqrt(pend(1,2)**2+pend(2,2)**2+pend(3,2)**2)
2705 pend(1,4)=xxm1pr(n,k)
2706 pend(2,4)=xym1pr(n,k)
2707 pend(3,4)=-xm1pr(n,k)*plc/2d0
2708 pend(4,4)=dsqrt(pend(1,4)**2+pend(2,4)**2+pend(3,4)**2)
2709 pend(1,3)=xxm2pr(n,k)
2710 pend(2,3)=xym2pr(n,k)
2711 pend(3,3)=-xm2pr(n,k)*plc/2d0
2712 pend(4,3)=dsqrt(pend(1,3)**2+pend(2,3)**2+pend(3,3)**2)
2716 pstg(1,1)=xxp1pr(n,k)+xxm2pr(n,k)
2717 pstg(2,1)=xyp1pr(n,k)+xym2pr(n,k)
2718 pstg(3,1)=(xp1pr(n,k)-xm2pr(n,k))*plc/2d0
2719 pstg(4,1)=(xp1pr(n,k)+xm2pr(n,k))*plc/2d0
2720 pstg(5,1)=dsqrt((pstg(4,1)-pstg(3,1))*(pstg(4,1)+pstg(3,1))
2721 & -pstg(1,1)**2-pstg(2,1)**2)
2722 pstg(1,2)=xxp2pr(n,k)+xxm1pr(n,k)
2723 pstg(2,2)=xyp2pr(n,k)+xym1pr(n,k)
2724 pstg(3,2)=(xp2pr(n,k)-xm1pr(n,k))*plc/2d0
2725 pstg(4,2)=(xp2pr(n,k)+xm1pr(n,k))*plc/2d0
2726 pstg(5,2)=dsqrt((pstg(4,2)-pstg(3,2))*(pstg(4,2)+pstg(3,2))
2727 & -pstg(2,2)**2-pstg(1,2)**2)
2733 if(ntry.gt.100)goto1001
2755 if(ish.ge.5)write(ifch,'(a,3x,2i7,i9)')' proj: '
2756 * ,(icp(l),l=1,2),idpj0
2757 if(ish.ge.5)write(ifch,'(a,3x,2i7,i9)')' targ: '
2758 * ,(ict(l),l=1,2),idtg0
2760 c determine string flavors
2762 call fstrfl(icp,ict,icp1,icp2,icm1,icm2
2763 * ,idp1pr(n,k),idp2pr(n,k),idm1pr(n,k),idm2pr(n,k)
2764 * ,iabs(idfpr(n,k)),iret)
2766 jerr(1)=jerr(1)+1 ! > 9 quarks per flavor attempted.
2767 ! OK when happens rarely.
2771 c check mass string 1
2773 ic(1)=icp1(1)+icm2(1)
2774 ic(2)=icp1(2)+icm2(2)
2775 if(ic(1).gt.0.or.ic(2).gt.0)then
2777 call iddeco(icp1,jcp1)
2778 call iddeco(icm2,jcm2)
2779 ammns=utamnx(jcp1,jcm2)
2780 if(ish.ge.7)write(ifch,'(a,2i7,2e12.3)')
2781 * ' string 1 - ic,mass,min.mass:',ic,am,ammns
2782 if(am.lt.ammns*facmss)then
2783 goto 777 !avoid virpom
2785 idend(1)=idtra(icp1,0,0,3)
2786 idend(3)=idtra(icm2,0,0,3)
2787 if(ish.ge.7)write(ifch,'(a,2i4)') ' string 1 - SE-ids:'
2788 * ,idend(1),idend(3)
2791 c check mass string 2
2793 ic(1)=icp2(1)+icm1(1)
2794 ic(2)=icp2(2)+icm1(2)
2795 if(ic(1).gt.0.or.ic(2).gt.0)then
2797 call iddeco(icp2,jcp2)
2798 call iddeco(icm1,jcm1)
2799 ammns=utamnx(jcp2,jcm1)
2800 if(ish.ge.7)write(ifch,'(a,2i7,2e12.3)')
2801 * ' string 2 - ic,mass,min.mass:',ic,am,ammns
2802 if(am.lt.ammns*facmss)then
2803 goto 777 !avoid virpom
2805 idend(2)=idtra(icp2,0,0,3)
2806 idend(4)=idtra(icm1,0,0,3)
2807 if(ish.ge.7)write(ifch,'(a,2i4)') ' string 2 - SE-ids:'
2808 * ,idend(2),idend(4)
2812 write(ifch,'(a,i10)')' pom: '
2814 write(ifch,'(a,2i5)')' str 1: '
2815 * ,idend(1),idend(3)
2816 write(ifch,'(a,2i5)')' str 2: '
2817 * ,idend(2),idend(4)
2818 write(ifch,'(a,2i7,1x,a)')' proj: '
2820 write(ifch,'(a,2i7,1x,a)')' targ: '
2831 call idtr4(idptl(ip),icini) !excited remnant ?
2832 if(ish.ge.5)write(ifch,*)'icini proj',icini
2833 & ,(icp(1)-icini(1)),(icp(2)-icini(2))
2834 if((icp(1)-icini(1))+(icp(2)-icini(2)).ne.0)iep(ip)=1
2835 call idtr4(idptl(maproj+it),icini)
2836 if(ish.ge.5)write(ifch,*)'icini targ',icini
2837 & ,(ict(1)-icini(1)),(ict(2)-icini(2))
2838 if((ict(1)-icini(1))+(ict(2)-icini(2)).ne.0)iet(it)=1
2839 if(ish.ge.5)write(ifch,*)'iep,iet ',iep(ip),iet(it)
2841 c write strings to /cptl/
2843 its=idp1pr(n,k)+idm2pr(n,k)
2844 call fstrwr(1,1,3,k,n)
2845 its=idp2pr(n,k)+idm1pr(n,k)
2846 call fstrwr(2,2,4,k,n)
2852 call utprix('ProSeF',ish,ishini,6)
2860 c-----------------------------------------------------------------------
2861 subroutine fstrfl(icp,ict,icp1,icp2,icm1,icm2
2862 * ,idp1,idp2,idm1,idm2,idfp,iret)
2863 c-----------------------------------------------------------------------
2864 c knowing the string end types (idp1,idp2,idm1,idm2)
2865 c and remnant flavors (icp,ict)
2866 c and remnant link of the string (idfp)
2867 c one determines quark flavors of string ends (icp1,icp2,icm1,icm2)
2868 c and updates remnant flavors (icp,ict)
2870 c iret=1 problem, more than 9 quarks per flavor attempted
2871 c-----------------------------------------------------------------------
2873 integer icp(2),ict(2),icp1(2),icp2(2),icm1(2),icm2(2)
2874 integer jcp(6,2),jct(6,2),jcpi(6,2),jcti(6,2)
2876 c data neuz/0/proz/0/dtaz/0/
2877 c save neuz,proz,dtaz
2879 call utpri('fstrfl',ish,ishini,7)
2890 if(idfp.eq.2)stop'fstrfl: should not happen (2). '
2891 if(idfp.eq.3)stop'fstrfl: should not happen (3). '
2892 if(idp1.eq.4)stop'fstrfl: diq code 4 not used any more'
2893 if(idm1.eq.4)stop'fstrfl: diq code 4 not used any more'
2894 if(idp2.eq.4)stop'fstrfl: diq code 4 not used any more'
2895 if(idm2.eq.4)stop'fstrfl: diq code 4 not used any more'
2896 if(idp1.eq.8)stop'fstrfl: fragm quarks not used any more'
2897 if(idp2.eq.8)stop'fstrfl: fragm quarks not used any more'
2898 if(idm1.eq.8)stop'fstrfl: fragm quarks not used any more'
2899 if(idm2.eq.8)stop'fstrfl: fragm quarks not used any more'
2901 c determine flavors of string ends (u,d,s)
2903 call iddeco(icp,jcpi)
2904 call iddeco(ict,jcti)
2905 call iddeco(icp,jcp)
2906 call iddeco(ict,jct)
2908 write(ifch,'(a,2i7,5x,6i2,3x,6i2,3x,i1)')' proj:',icp,jcp
2909 write(ifch,'(a,2i7,5x,6i2,3x,6i2,3x,i1)')' targ:',ict,jct
2934 iq(1,1)=idrafl(iclpro,jcp,1,'s',iret)
2938 iq(1,2)=idrafl(iclpro,jcp,2,'s',iret)
2942 iq(1,4)=idrafl(icltar,jct,1,'s',iret)
2946 iq(1,3)=idrafl(icltar,jct,2,'s',iret)
2953 iq(1,1)=idrafl(iclpro,jcp,1,'s',iret1)
2957 iq(1,4)=idrafl(icltar,jct,1,'s',iret4)
2961 iq(1,2)=idrafl(iclpro,jcp,2,'s',iret2)
2965 iq(1,3)=idrafl(icltar,jct,2,'s',iret3)
2969 c diquarks, code 5 (former valence, but actually sea)
2973 c iq(1,1)=idraflx(fc,iclpro,jcp,2,'s',iret)
2974 c if(iq(1,1).eq.3)fc=fc*puds
2975 c iq(2,1)=idraflx(fc,iclpro,jcp,2,'s',iret)
2976 iq(1,1)=idrafl(iclpro,jcp,2,'d',iret)
2977 iq(2,1)=idrafl(iclpro,jcp,2,'d',iret)
2981 c iq(1,4)=idraflx(fc,icltar,jct,2,'s',iret)
2982 c if(iq(1,4).eq.3)fc=fc*puds
2983 c iq(2,4)=idraflx(fc,icltar,jct,2,'s',iret)
2984 iq(1,4)=idrafl(icltar,jct,2,'d',iret)
2985 iq(2,4)=idrafl(icltar,jct,2,'d',iret)
2989 c iq(1,2)=idraflx(fc,iclpro,jcp,1,'s',iret)
2990 c if(iq(1,2).eq.3)fc=fc*puds
2991 c iq(2,2)=idraflx(fc,iclpro,jcp,1,'s',iret)
2992 iq(1,2)=idrafl(iclpro,jcp,1,'d',iret)
2993 iq(2,2)=idrafl(iclpro,jcp,1,'d',iret)
2997 c iq(1,3)=idraflx(fc,icltar,jct,1,'s',iret)
2998 c if(iq(1,3).eq.3)fc=fc*puds
2999 c iq(2,3)=idraflx(fc,icltar,jct,1,'s',iret)
3000 iq(1,3)=idrafl(icltar,jct,1,'d',iret)
3001 iq(2,3)=idrafl(icltar,jct,1,'d',iret)
3004 if(iret.ne.0)goto 1000
3007 c in case of saturated remnants, use the same flavor for quark and anti-quark
3009 if(iret1.ne.0.and.iret2.ne.0)then
3010 call iddeco(icp,jcp)
3011 if(rangen().gt.0.5)then
3016 elseif(iret1.eq.0.and.iret2.ne.0.and.idp1.eq.1)then
3017 call iddeco(icp,jcp)
3019 if(idp1.eq.4)iq(2,1)=iq(1,1)
3020 c if(idp2.eq.4)iq(2,2)=iq(1,2)
3021 elseif(iret2.eq.0.and.iret1.ne.0.and.idp2.eq.1)then
3022 call iddeco(icp,jcp)
3024 if(idp1.eq.4)iq(2,1)=iq(1,1)
3025 c if(idp2.eq.4)iq(2,2)=iq(1,2)
3026 elseif(iret1.ne.0.or.iret2.ne.0)then
3031 if(iret3.ne.0.and.iret4.ne.0)then
3032 call iddeco(ict,jct)
3033 if(rangen().gt.0.5)then
3038 elseif(iret3.eq.0.and.iret4.ne.0.and.idm1.eq.1)then
3039 call iddeco(ict,jct)
3041 c if(idm2.eq.4)iq(2,3)=iq(1,3)
3042 if(idm1.eq.4)iq(2,4)=iq(1,4)
3043 elseif(iret4.eq.0.and.iret3.ne.0.and.idm2.eq.1)then
3044 call iddeco(ict,jct)
3046 c if(idm2.eq.4)iq(2,3)=iq(1,3)
3047 if(idm1.eq.4)iq(2,4)=iq(1,4)
3048 elseif(iret3.ne.0.or.iret4.ne.0)then
3056 call idenco(jcp,icp,iret)
3057 if(iret.ne.0)goto 1000
3058 call idenco(jct,ict,iret)
3059 if(iret.ne.0)goto 1000
3065 if(ish.ge.7)write(ifch,'(a,2i3,4x,2i3)')
3066 *' string 1, string ends:',ifla,iflb,iflc,ifld
3070 icp1(1)=10**(6-ifla)
3074 icp1(2)=10**(6-ifla)
3075 icp1(2)=icp1(2)+10**(6-iflb)
3082 icm2(2)=10**(6-iflc)
3084 icm2(1)=10**(6-iflc)
3085 icm2(1)=icm2(1)+10**(6-ifld)
3094 if(ish.ge.7)write(ifch,'(a,2i3,4x,2i3)')
3095 *' string 2, string ends:',ifla,iflb,iflc,ifld
3099 icm1(1)=10**(6-ifla)
3103 icm1(2)=10**(6-ifla)
3104 icm1(2)=icm1(2)+10**(6-iflb)
3111 icp2(2)=10**(6-iflc)
3113 icp2(1)=10**(6-iflc)
3114 icp2(1)=icp2(1)+10**(6-ifld)
3120 write(ifch,'(a,2i7,4x,2i7)')
3121 * ' SE-forw:',icp1(1),icp1(2),icp2(1),icp2(2)
3122 write(ifch,'(a,2i7,4x,2i7)')
3123 * ' SE-back:',icm1(1),icm1(2),icm2(1),icm2(2)
3124 write(ifch,'(a,2i7,5x,6i2,3x,6i2)')' proj:',icp,jcp
3125 write(ifch,'(a,2i7,5x,6i2,3x,6i2)')' targ:',ict,jct
3132 call utprix('fstrfl',ish,ishini,7)
3136 c-----------------------------------------------------------------------
3137 integer function jdrafl(icl,jc,mod,iret)
3138 c-----------------------------------------------------------------------
3140 c returns random flavor of a quark
3143 c jc : quark content of remnant
3144 c returns random flavor and update remant with corresponding q-qbar pair \
3145 c if there is enough place (else iret=1)
3147 c id=1 u, id=2 d, id=3 s
3148 c-----------------------------------------------------------------------
3152 c write(*,*)'entry jdrafl, j,c,jc: ',j,c,jc
3161 if(r.gt.(pu+pd).and.ps.gt.1d-10)then
3163 elseif(r.gt.pu.and.pd.gt.1d-10)then
3169 i=1+int((2.+rstras(icl))*rangen())
3173 c write(*,*)'jc before updating',jc
3174 c write(*,*)'i,j,jc',i,j,jc
3177 call idsufl2(i,1,jc,iret)
3178 call idsufl2(i,2,jc,iret)
3185 cc-----------------------------------------------------------------------
3186 c subroutine fremfl(icp,ict,iret)
3187 cc-----------------------------------------------------------------------
3188 cc checks projectile and target flavor (icp,ict)
3189 cc in case of reggeon exchange they do not correspond to hadrons.
3190 cc one transfers therefore flavor from one side to the other in order
3191 cc to have hadron flavor.
3192 cc icp and ict are modified correspondingly
3193 cc-----------------------------------------------------------------------
3194 c include 'epos.inc'
3195 c integer icp(2),ict(2),jcp(6,2),jct(6,2),kp(4),kt(4)
3197 c call utpri('fremfl',ish,ishini,7)
3204 c call iddeco(icp,jcp)
3205 c call iddeco(ict,jct)
3212 c kp(l)=jcp(l,1)-jcp(l,2)
3213 c kt(l)=jct(l,1)-jct(l,2)
3214 c iakp=iakp+iabs(kp(l))
3215 c iakt=iakt+iabs(kt(l))
3219 c if(ish.ge.7)write(ifch,*)'iak_p:',iakp,' ik_p:',ikp
3220 c if(ish.ge.7)write(ifch,*)'iak_t:',iakt,' ik_t:',ikt
3223 c if(ikp.eq.4.or.ikp.eq.-2)then
3224 c ifl=idrafl(jcp,1,'v',iret)
3225 c iqp=2 ! subtract quark
3227 c elseif(ikp.eq.-4.or.ikp.eq.2)then
3228 c ifl=idrafl(jcp,2,'v',iret)
3229 c iqp=1 ! subtract antiquark
3230 c iqt=2 ! add antiquark
3232 c call utstop('fremfl&')
3234 c elseif(iakt.eq.4)then
3235 c if(ikt.eq.4.or.ikt.eq.-2)then
3236 c ifl=idrafl(jct,1,'v',iret)
3238 c iqt=2 ! subtract quark
3239 c elseif(ikt.eq.-4.or.ikt.eq.2)then
3240 c ifl=idrafl(jct,2,'v',iret)
3241 c iqp=2 ! add antiquark
3242 c iqt=1 ! subtract antiquark
3244 c call utstop('fremfl&')
3246 c elseif(iakp.eq.3)then
3248 c ifl=idrafl(jcp,1,'v',iret)
3249 c iqp=2 ! subtract quark
3252 c ifl=idrafl(jcp,2,'v',iret)
3253 c iqp=1 ! subtract antiquark
3254 c iqt=2 ! add antiquark
3256 c elseif(iakt.eq.3)then
3258 c ifl=idrafl(jct,1,'v',iret)
3260 c iqt=2 ! subtract quark
3262 c ifl=idrafl(jct,2,'v',iret)
3263 c iqp=2 ! add antiquark
3264 c iqt=1 ! subtract antiquark
3266 c elseif(iakp.eq.2)then
3268 c ifl=idrafl(jct,1,'v',iret)
3270 c iqt=2 ! subtract quark
3272 c ifl=idrafl(jct,2,'v',iret)
3273 c iqp=2 ! add antiquark
3274 c iqt=1 ! subtract antiquark
3276 c elseif(iakt.eq.2)then
3278 c ifl=idrafl(jct,1,'v',iret)
3279 c iqp=2 ! subtract quark
3282 c ifl=idrafl(jct,2,'v',iret)
3283 c iqp=1 ! subtract antiquark
3284 c iqt=2 ! add antiquark
3286 c elseif(iakp.eq.1)then
3288 c ifl=idrafl(jcp,2,'v',iret)
3289 c iqp=2 ! add antiquark
3290 c iqt=1 ! subtract antiquark
3292 c ifl=idrafl(jcp,1,'v',iret)
3294 c iqt=2 ! subtract quark
3296 c elseif(iakt.eq.1)then
3298 c ifl=idrafl(jct,2,'v',iret)
3299 c iqp=1 ! subtract antiquark
3300 c iqt=2 ! add antiquark
3302 c ifl=idrafl(jct,1,'v',iret)
3303 c iqp=2 ! subtract quark
3307 c call utstop('fremfl: error&')
3310 c if(ish.ge.7)write(ifch,*)'iq_p:',iqp,' iq_t:',iqt,' if:',ifl
3311 c call uticpl(icp,ifl,iqp,iret)
3312 c if(iret.ne.0)goto1000
3313 c call uticpl(ict,ifl,iqt,iret)
3314 c if(iret.ne.0)goto1000
3320 c call utprix('fremfl',ish,ishini,7)
3324 c-----------------------------------------------------------------------
3325 subroutine fstrwr(j,ii,jj,k,n)
3326 c-----------------------------------------------------------------------
3327 c take pstg(5,j),pend(4,ii),idend(ii),pend(4,jj),idend(jj) (/cems/)
3328 c and write it to /cptl/
3329 c-----------------------------------------------------------------------
3331 c ii,jj: string end (1,2: proj; 3,4: targ)
3332 c k: current collision
3333 c n: current pomeron
3334 c-----------------------------------------------------------------------
3337 include 'epos.incems'
3339 double precision pstg,pend
3340 common/cems/pstg(5,2),pend(4,4),idend(4)
3341 common/emsptl/nppr(npommx,kollmx),npproj(mamx),nptarg(mamx)
3342 double precision pp(4)
3345 call utpri('fstrwr',ish,ishini,7)
3347 if(idend(ii).ne.0.and.idend(jj).ne.0)then
3351 call utlob2(1,pstg(1,j),pstg(2,j),pstg(3,j),pstg(4,j),pstg(5,j)
3352 * ,pend(1,ii),pend(2,ii),pend(3,ii),pend(4,ii),20)
3355 pp(3)=.5d0*pstg(5,j)
3356 pp(4)=.5d0*pstg(5,j)
3358 * (-1,pend(1,ii),pend(2,ii),pend(3,ii),pp(1),pp(2),pp(3))
3359 call utlob2(-1,pstg(1,j),pstg(2,j),pstg(3,j),pstg(4,j),pstg(5,j)
3360 * ,pp(1),pp(2),pp(3),pp(4),21)
3363 if(ifrptl(1,npom).eq.0)ifrptl(1,npom)=nptl+1
3364 ifrptl(2,npom)=nptl+2
3368 pptl(1,nptl)=sngl(pp(1))
3369 pptl(2,nptl)=sngl(pp(2))
3370 pptl(3,nptl)=sngl(pp(3))
3371 pptl(4,nptl)=sngl(pp(4))
3378 xorptl(1,nptl)=coord(1,k)
3379 xorptl(2,nptl)=coord(2,k)
3380 xorptl(3,nptl)=coord(3,k)
3381 xorptl(4,nptl)=coord(4,k)
3382 tivptl(1,nptl)=xorptl(4,nptl)
3383 tivptl(2,nptl)=xorptl(4,nptl)
3384 idptl(nptl)=idend(ii)
3385 ityptl(nptl)=ityptl(npom)+j
3393 pptl(i,nptl)=sngl(pstg(i,j))-pptl(i,nptl-1)
3398 iorptl(nptl)=nppr(n,k)
3402 xorptl(1,nptl)=coord(1,k)
3403 xorptl(2,nptl)=coord(2,k)
3404 xorptl(3,nptl)=coord(3,k)
3405 xorptl(4,nptl)=coord(4,k)
3406 tivptl(1,nptl)=xorptl(4,nptl)
3407 tivptl(2,nptl)=xorptl(4,nptl)
3408 idptl(nptl)=idend(jj)
3409 ityptl(nptl)=ityptl(npom)+j
3416 write(ifch,100)' kink:',(pptl(l,nptl-1),l=1,4),idptl(nptl-1)
3417 write(ifch,100)' kink:',(pptl(l,nptl),l=1,4),idptl(nptl)
3420 elseif(idend(ii).ne.0.and.idend(jj).eq.0)then
3425 if(ifrptl(1,npom).eq.0)ifrptl(1,npom)=nptl+1
3426 ifrptl(2,npom)=nptl+1
3430 idptl(nptl)=idend(ii)
3431 pptl(1,nptl)=sngl(pstg(1,j))
3432 pptl(2,nptl)=sngl(pstg(2,j))
3433 pptl(3,nptl)=sngl(pstg(3,j))
3434 pptl(4,nptl)=sngl(pstg(4,j))
3435 pptl(5,nptl)=sngl(pstg(5,j))
3441 xorptl(1,nptl)=coord(1,k)
3442 xorptl(2,nptl)=coord(2,k)
3443 xorptl(3,nptl)=coord(3,k)
3444 xorptl(4,nptl)=coord(4,k)
3445 tivptl(1,nptl)=coord(4,k)
3446 call idtau(idptl(nptl),pptl(4,nptl),pptl(5,nptl),taugm)
3447 tivptl(2,nptl)=tivptl(1,nptl)+taugm*(-alog(rangen()))
3448 ityptl(nptl)=ityptl(npom)+2+j
3455 write(ifch,100)' res:',(pptl(l,nptl),l=1,4),idptl(nptl)
3457 elseif(idend(ii).eq.0.and.idend(jj).eq.0)then
3460 call utstop('error in fstrwr&')
3463 100 format(a,4e9.3,i5)
3466 call utprix('fstrwr',ish,ishini,7)
3470 c-----------------------------------------------------------------------
3471 subroutine ProReF(ir,m)
3472 c-----------------------------------------------------------------------
3473 c proposes flavor for remnant m for proj (ir=1) or target (ir=-1)
3474 c and writes remnant into /cptl/ as string or hadron
3475 c ityptl definitions:
3476 c 51 41 ... rmn drop
3477 c 52 42 ... rmn str inel
3478 c 53 43 ... rmn str diff
3479 c 54 44 ... rmn str after droplet or hadron split
3481 c 56 46 ... rmn res after droplet or hadron split
3482 c 57 47 ... rmn res after all Pomeron killed
3483 c 58 48 ... rmn res from diff
3484 c 59 49 ... hadron split
3485 c-----------------------------------------------------------------------
3488 include 'epos.incems'
3490 double precision plc,s !,ptt1,ptt2
3492 common/cdfptl/idfptl(mxptl)
3493 double precision tpro,zpro,ttar,ztar,ttaus,detap,detat,zor,tor
3494 common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat
3495 common /cncl/xproj(mamx),yproj(mamx),zproj(mamx)
3496 * ,xtarg(mamx),ytarg(mamx),ztarg(mamx)
3497 double precision amasmin,amasini,mdrmax
3498 integer icf(2),icb(2)
3499 integer jcf(nflav,2),jcdummy(nflav,2)
3500 logical gdrop, ghadr,gproj,gtarg
3501 double precision ept(5),ep(4),aa(5),am2t
3502 common/emsptl/nppr(npommx,kollmx),npproj(mamx),nptarg(mamx)
3503 common /ems12/iodiba,bidiba ! defaut iodiba=0. if iodiba=1, study H-Dibaryon
3505 call utpri('ProReF',ish,ishini,3)
3507 if(ir.ne.1.and.ir.ne.-1)stop'ProReF: wrong ir'
3517 c if(kolp(m).le.0)goto1000
3518 if(iep(m).le.-1)goto1000
3525 elseif(ir.eq.-1)then
3526 c if(kolt(m).le.0)goto1000
3527 if(iet(m).le.-1)goto1000
3535 call utstop('ProReF: ir ???&')
3537 if(ish.ge.3)write(ifch,*)'remnant particle index:',mm
3539 if(ish.ge.8)call alist('ProRef&',1,nptl)
3544 minfra=min(minfra,nptlini) !for trigger condition
3547 ept(l)=dble(pptl(l,mm))
3553 c initialize forward and backward ic (to transform remnant into string)
3565 call iddeco(icf,jcf)
3566 call idquacjc(jcf,nqu,naq)
3570 amasmin=dble(fremnux2(icf))**2.d0
3571 if(ept(5).le.0.)then
3572 ept(5)=sqrt(2*amasmin)
3574 call utmsg('ProReF')
3575 write(ifch,*)'zero remnant mass -> amasmin'
3579 am2t=(ept(4)+ept(3))*(ept(4)-ept(3))-(ept(1)**2+ept(2)**2)
3581 & .and.(am2t.le.0d0.or.abs(am2t-ept(5)*ept(5)).gt.ept(5)))then
3582 write(ifch,*)'Precision problem in ProRef, p:',
3583 & (ept(k),k=1,4),ept(5)*ept(5),am2t
3585 ept(4)=sqrt(ept(3)*ept(3)+ept(2)*ept(2)+ept(1)*ept(1)
3590 * write(ifch,'(a,5e11.3,2i7)')' proj:'
3591 * ,(sngl(ept(k)) ,k=1,5),(icproj(k,m) ,k=1,2)
3593 * write(ifch,'(a,5e11.3,2i7)')' targ:'
3594 * ,(sngl(ept(k)) ,k=1,5),(ictarg(k,m),k=1,2)
3597 amasini=ept(5)*ept(5)
3599 c mdrmax=amasmin+dble(amdrmax*amdrmax)
3600 mdrmax=dble(fremnux(icf)+amdrmax)**2.d0
3602 if(ish.ge.4)write(ifch,*)'remnant masses:',am2t,amasini,amasmin
3605 c.............................exotic ...................................
3607 c if(amasini.gt.amasmin.and.irmdropx.eq.1)then
3610 c & .not.((nqu.eq.3.and.naq.eq.0).or.(nqu.eq.0.and.naq.eq.3)
3611 if(.not.((nqu.eq.3.and.naq.eq.0).or.(nqu.eq.0.and.naq.eq.3)
3612 & .or.(nqu.eq.1.and.naq.eq.1))
3613 & .and.amasini.gt.amasmin.and.irmdropx.eq.1)then
3616 c & .not.((nqu.eq.3.and.naq.eq.0).or.(nqu.eq.0.and.naq.eq.3)
3617 c & .or.(nqu.eq.1.and.naq.eq.1)).or.
3618 c & (iept.ne.0.and.iept.le.2.and.reminv/ept(5).gt.rangen()))
3619 c & .and.amasini.gt.amasmin.and.irmdropx.eq.1)then
3621 !print*,'-------------------------------------------' !!!
3623 !print*,icf,sqrt(amasini),sqrt(amasmin),sqrt(mdrmax) !!!
3625 if(amasini.gt.mdrmax.or.(jcf(4,1)+jcf(4,2).ne.0))then !charm not possible in droplet
3626 call getdroplet(ir,icf,jcf,ept,aa,gdrop,mdrmax)
3627 !--------------------------------
3628 !emit a droplet, update the remnant string flavour and 5-momentum
3630 ! ir ......... 1 projectile, -1 target remnant
3631 ! ept ........ remnant 5-momentum
3632 ! jcf ........ remnant jc
3634 ! gdrop ... .true. = successful droplet emission
3635 ! icf, ept ....... droplet ic and 5-momentum
3636 ! jcf, a ......... remnant string jc and 5-momentum
3637 ! .false. = unsuccessful
3638 ! jcf, ept .... unchanged,
3639 ! emits hadrons instead of droplet
3640 c ! considered as droplet jc and 5-momentum
3641 !-------------------------------------
3642 if(.not.gdrop)goto 500
3646 !also in case of unsuccessful drop emission, then remnant = droplet !
3654 c Remnant radius to have eps=dens GeV/fm3
3655 radptl(nptl)=(3.*sngl(ept(5))/4./pi/dens)**0.3333
3658 pptl(l,nptl)=sngl(ept(l))
3660 idx=idtra(icf,0,0,3)
3663 call idres(idx,amx,idrx,iadjx,1)
3668 idptl(nptl)=8*10**8+icf(1)*100+icf(2)/100
3680 if(iept.eq.6)ityptl(nptl)=47
3683 if(iept.eq.6)ityptl(nptl)=57
3690 xorptl(1,nptl)=xorptl(1,mm)
3691 xorptl(2,nptl)=xorptl(2,mm)
3692 xorptl(3,nptl)=xorptl(3,mm)
3695 call idtau(idptl(nptl),pptl(4,nptl),pptl(5,nptl),taugm)
3696 tivptl(2,nptl)=tivptl(1,nptl)+taugm*(-alog(rangen()))
3701 if(ish.ge.3)write(ifch,*)'Proref,ept(5),id',ept(5),idptl(nptl)
3702 !print*,nptl,idptl(nptl),sngl(ept(5)),pptl(5,nptl) !!!
3704 !..........remnant update
3705 if(gdrop)then !drop emission: new remnant -> ept, icf
3710 call idquacjc(jcf,nqu,naq)
3711 call idenco(jcf,icf,iret)
3712 if(iret.eq.1)call utstop('Pb in ProRef in strg+drop process&')
3713 !!! print*,'new remnant:',icf,ept(5) !!!
3718 pptl(l,nptl)=sngl(ept(l))
3720 idptl(nptl)=idptl(mm)
3726 xorptl(1,nptl)=xorptl(1,mm)
3727 xorptl(2,nptl)=xorptl(2,mm)
3728 xorptl(3,nptl)=xorptl(3,mm)
3731 tivptl(2,nptl)=ainfin
3742 !........decay mini-droplet......
3745 if(iabs(idptl(mm)).gt.10**8)then
3746 if(ish.ge.3)write(ifch,*)'Make droplet'
3747 if(nptlb.gt.mxptl-10)call utstop('ProRef: mxptl too small&')
3749 if(ifrade.gt.0.and.ispherio.eq.0)call hnbaaa(mm,iret)!Decay remn
3750 if(iret.ne.1.and.nptl.ne.nptlb)then ! ---successful decay---
3751 istptl(mm)=istptl(mm)+1
3752 ifrptl(1,mm)=nptlb+1
3755 x=xorptl(1,mm)+(t-xorptl(4,mm))*pptl(1,mm)/pptl(4,mm)
3756 y=xorptl(2,mm)+(t-xorptl(4,mm))*pptl(2,mm)/pptl(4,mm)
3757 z=xorptl(3,mm)+(t-xorptl(4,mm))*pptl(3,mm)/pptl(4,mm)
3758 do 21 n=nptlb+1,nptl
3764 if(idfptl(mm).eq.0)then
3769 radius=0.8*sqrt(rangen())
3773 xorptl(1,n)=x + radius*cos(phi)
3774 xorptl(2,n)=y + radius*sin(phi)
3778 zor=dble(xorptl(3,iioo))
3779 tor=dble(xorptl(4,iioo))
3780 call idquac(iioo,nq,ndummy1,ndummy2,jcdummy)
3782 tauran=-taurea*alog(r)
3783 call jtaix(n,tauran,zor,tor,zis,tis)
3784 tivptl(1,n)=amax1(ti,tis)
3785 call idtau(idptl(n),pptl(4,n),pptl(5,n),taugm)
3787 tivptl(2,n)=t+taugm*(-alog(r))
3790 if(iept.eq.6)ityptl(n)=47
3793 if(iept.eq.6)ityptl(n)=57
3800 if(iabs(idptl(nptlb+1)).le.6) then
3802 if(ish.ge.1)write (ifmt,*)'string from drop:nptlb+1,nptl:'
3813 elseif(ifrade.gt.0.and.ispherio.eq.0)then ! Unsuccessful decay
3815 if(ish.ge.4)write(ifch,*)
3816 * '***** Unsuccessful remnant cluster decay'
3817 * ,' --> do RemoveHadrons instead.'
3825 if(idrop.eq.1)goto 1000
3826 !successful drop decay, no additional string, nothing to do
3830 c...............................................................
3833 if(gdrop)mm=nptlini+2
3837 c........................remove hadrons.........................
3842 if(.not.((nqu.eq.3.and.naq.eq.0).or.(nqu.eq.0.and.naq.eq.3)
3843 & .or.(nqu.eq.1.and.naq.eq.1)))then
3844 if(irmdropx.eq.irmdrop)then
3846 !call utmsg('ProReF')
3847 !write(ifch,*)'***** condition for droplet treatment: '
3848 !write(ifch,*)'***** amasini.gt.amasmin.and.irmdropx.eq.1 = '
3849 !* ,amasini.gt.amasmin.and.irmdropx.eq.1
3850 !write(ifch,*)'***** amasini,amasmin,irmdropx:'
3851 !* ,amasini,amasmin,irmdropx
3852 !write(ifch,*)'***** nqu,naq:',nqu,naq
3853 !write(ifch,*)'***** call RemoveHadrons'
3856 call RemoveHadrons(gproj,gtarg,ghadr,m,mm,jcf,icf,ept)
3859 c........................ determine idr (0=string, else=resonance).......
3861 if(icf(1).eq.0.and.icf(2).eq.0)then
3868 call idres(id,am,idr,iadj,1)
3869 if(iabs(mod(idr,10)).le.2.and.idr.ne.0)then
3873 endif !ckeck on-shell mass (see uti)
3874 if(iadj.ne.0.and.iept.gt.0.and.ept(5).gt.0.d0
3875 & .and.(dabs((ept(4)+ept(3))*(ept(4)-ept(3))
3876 $ -ept(2)**2-ept(1)**2-dble(am)**2).gt.0.3d0))idr=0
3879 write(ifch,'(a,5e11.3)')' updt:',(sngl(ept(k)) ,k=1,5)
3880 write(ifch,*)' icf: ',icf,' idr: ',idr,' iept: ',iept
3883 if(iept.eq.3)stop'ProReF: iept=3 ???'
3885 c...........................................string...................
3886 if(iept.gt.0.and.idr.eq.0)then
3888 !... nqu of remainder string
3891 if(gdrop)anstrg1=anstrg1+1
3893 call iddeco(icf,jcf)
3896 nqu=nqu+jcf(l,1)-jcf(l,2)
3899 if(zbarfl.lt.0.)stop'ProReF: not supported any more. '
3901 !......determine forward momentum ep
3903 !ptt=0.5*min(zopmax,zopinc*zz)
3905 !ptt1=dble(ptt*cos(phi))
3906 !ptt2=dble(ptt*sin(phi))
3910 ep(3)=ir*0.5d0*ept(5)
3913 call utlob2(-1,ept(1),ept(2),ept(3),ept(4),ept(5)
3914 * ,ep(1),ep(2),ep(3),ep(4),25)
3916 !....determine forward and backward flavor icf, icb
3919 c if(iept.le.2.and.ept(5)/reminv.lt.rangen())ireminv=1
3920 c if(iept.eq.2.and.ept(5).lt.reminv.and.rangen().lt.0.5)ireminv=1
3921 if(iept.eq.6.and.rangen().lt.0.25)ireminv=1
3923 c if(reminv/ept(5).gt.rangen())ireminv=1
3924 c elseif(iept.eq.6)then
3927 if(nqu.eq.3)then !---baryon---
3928 iq=idrafl(iclpt,jcf,1,'v',iret)
3929 call uticpl(icf,iq,2,iret) ! antiquark
3930 call uticpl(icb,iq,1,iret) ! quark
3931 if(ireminv.eq.1)then
3932 iq=idrafl(iclpt,jcf,1,'v',iret)
3933 call uticpl(icf,iq,2,iret) ! antiquark
3934 call uticpl(icb,iq,1,iret) ! quark
3936 elseif(nqu.eq.-3)then !---antibaryon---
3937 iq=idrafl(iclpt,jcf,2,'v',iret)
3938 call uticpl(icf,iq,1,iret) ! quark
3939 call uticpl(icb,iq,2,iret) ! antiquark
3940 if(ireminv.eq.1)then
3941 iq=idrafl(iclpt,jcf,2,'v',iret)
3942 call uticpl(icf,iq,1,iret) ! quark
3943 call uticpl(icb,iq,2,iret) ! antiquark
3945 elseif(nqu.eq.0)then !---meson---
3946 iq1=idrafl(iclpt,jcf,1,'v',iret)
3947 iq2=idrafl(iclpt,jcf,2,'v',iret)
3948 if(rangen().gt.0.5)then
3949 call uticpl(icf,iq1,2,iret) ! subtract quark
3950 call uticpl(icb,iq1,1,iret) ! add quark
3952 call uticpl(icf,iq2,1,iret) ! subtract antiquark
3953 call uticpl(icb,iq2,2,iret) ! add antiquark
3955 c elseif(nqu.eq.0)then !---meson---
3956 c if(iept.ne.1.and.iept.ne.6.and.rangen().lt.0.5)then
3957 c iq=idrafl(iclpt,jcf,1,'v',iret)
3958 c call uticpl(icf,iq,2,iret) ! subtract quark
3959 c call uticpl(icb,iq,1,iret) ! add quark
3961 cc put quark in forward direction always for inelastic
3962 c iq=idrafl(iclpt,jcf,2,'v',iret)
3963 c call uticpl(icf,iq,1,iret) ! subtract antiquark
3964 c call uticpl(icb,iq,2,iret) ! add antiquark
3967 call utmsg('ProReF')
3968 write(ifch,*)'***** neither baryon nor antibaryon nor meson.'
3969 write(ifch,*)'***** number of net quarks:',nqu
3970 call utstop('ProRef&')
3972 c if(nqu.eq.3)then !---baryon---
3973 c iq1=idrafl(iclpt,jcf,1,'v',iret)
3974 c iq2=idrafl(iclpt,jcf,1,'v',iret)
3975 c iq3=idrafl(iclpt,jcf,1,'v',iret)
3976 c amdqa=qmass(iq2)+qmass(iq3)+qmass(0)
3977 c if(rangen().lt.qmass(iq1)/amdqa)ireminv=1
3978 c if(ireminv.ne.1)then
3979 c call uticpl(icf,iq1,2,iret) ! antiquark
3980 c call uticpl(icb,iq1,1,iret) ! quark
3982 c call uticpl(icf,iq2,2,iret) ! antiquark
3983 c call uticpl(icb,iq2,1,iret) ! quark
3984 c call uticpl(icf,iq3,2,iret) ! antiquark
3985 c call uticpl(icb,iq3,1,iret) ! quark
3987 c elseif(nqu.eq.-3)then !---antibaryon---
3988 c iq1=idrafl(iclpt,jcf,2,'v',iret)
3989 c iq2=idrafl(iclpt,jcf,2,'v',iret)
3990 c iq3=idrafl(iclpt,jcf,2,'v',iret)
3991 c amdqa=qmass(iq2)+qmass(iq3)+qmass(0)
3992 c if(rangen().lt.qmass(iq1)/amdqa)ireminv=1
3993 c if(ireminv.ne.1)then
3994 c call uticpl(icf,iq1,1,iret) ! antiquark
3995 c call uticpl(icb,iq1,2,iret) ! quark
3997 c call uticpl(icf,iq2,1,iret) ! antiquark
3998 c call uticpl(icb,iq2,2,iret) ! quark
3999 c call uticpl(icf,iq3,1,iret) ! antiquark
4000 c call uticpl(icb,iq3,2,iret) ! quark
4002 c elseif(nqu.eq.0)then !---meson---
4003 c iq1=idrafl(iclpt,jcf,1,'v',iret)
4004 c iq2=idrafl(iclpt,jcf,2,'v',iret)
4005 c if(rangen().lt.qmass(iq1)/qmass(iq2))then
4006 cc if(rangen().gt.0.5)then
4007 c call uticpl(icf,iq1,2,iret) ! subtract quark
4008 c call uticpl(icb,iq1,1,iret) ! add quark
4010 c call uticpl(icf,iq2,1,iret) ! subtract antiquark
4011 c call uticpl(icb,iq2,2,iret) ! add antiquark
4014 c call utmsg('ProReF')
4015 c write(ifch,*)'***** neither baryon nor antibaryon nor meson.'
4016 c write(ifch,*)'***** number of net quarks:',nqu
4017 c call utstop('ProRef&')
4020 !..... forward string end
4023 if(nptl.gt.mxptl)call utstop('ProRef: mxptl too small&')
4024 pptl(1,nptl)=sngl(ep(1))
4025 pptl(2,nptl)=sngl(ep(2))
4026 pptl(3,nptl)=sngl(ep(3))
4027 pptl(4,nptl)=sngl(ep(4))
4031 if(.not.gdrop)istptl(mm)=41
4033 if(nmes.eq.0.and.nbar.eq.0.and..not.gdrop)ifrptl(1,mm)=nptl
4035 xorptl(1,nptl)=xorptl(1,mm)
4036 xorptl(2,nptl)=xorptl(2,mm)
4037 xorptl(3,nptl)=xorptl(3,mm)
4038 xorptl(4,nptl)=xorptl(4,mm)
4039 tivptl(1,nptl)=xorptl(4,nptl)
4040 tivptl(2,nptl)=xorptl(4,nptl)
4041 idptl(nptl)=idtra(icf,0,0,3)
4043 if(iep(m).lt.1)stop'ProReF: iep(m)<1 '
4044 ityptl(nptl)=41+iep(m) ! =42 =43 =47
4045 if(gdrop.and.iep(m).ne.6)ityptl(nptl)=44
4046 if(ghadr)ityptl(nptl)=44
4048 if(iet(m).lt.1)stop'ProReF: iet(m)<1 '
4049 ityptl(nptl)=51+iet(m) !=52 =53 =57
4050 if(gdrop.and.iet(m).ne.6)ityptl(nptl)=54
4051 if(ghadr)ityptl(nptl)=54
4056 !write(6,'(a,i9,$)')' ',idptl(nptl) !======================
4063 write(ifch,'(a,5e11.3,$)')' kink:',(pptl(k,nptl),k=1,5)
4064 write(ifch,*)' id: ',idptl(nptl)
4066 !....... backward string end
4069 if(nptl.gt.mxptl)call utstop('ProRef: mxptl too small&')
4072 pptl(i,nptl)=sngl(ept(i)-ep(i))
4073 pptl2=pptl2+pptl(i,nptl)*pptl(i,nptl)
4075 pptl(4,nptl)=sqrt(pptl2)
4076 pptl2=sngl(ept(4)-ep(4))
4077 if(ish.ge.1.and.abs(pptl2-pptl(4,nptl)).gt.max(0.1,
4078 & 0.1*abs(pptl2)))then
4080 & 'Warning in ProRef: inconsistent backward string end energy !'
4081 & ,pptl(4,nptl),pptl2,abs(pptl2-pptl(4,nptl))
4082 if(ish.ge.2)write(ifch,*)
4083 & 'Warning in ProRef: inconsistent backward string end energy !'
4084 & ,(pptl(kkk,nptl),kkk=1,4),pptl2,abs(pptl2-pptl(4,nptl))
4093 xorptl(1,nptl)=xorptl(1,mm)
4094 xorptl(2,nptl)=xorptl(2,mm)
4095 xorptl(3,nptl)=xorptl(3,mm)
4096 xorptl(4,nptl)=xorptl(4,mm)
4097 tivptl(1,nptl)=xorptl(4,nptl)
4098 tivptl(2,nptl)=xorptl(4,nptl)
4099 idptl(nptl)=idtra(icb,0,0,3)
4101 ityptl(nptl)=41+iep(m) ! =42 =43 =47
4102 if(gdrop.and.iep(m).ne.6)ityptl(nptl)=44
4103 if(ghadr)ityptl(nptl)=44
4105 ityptl(nptl)=51+iet(m) !=52 =53 =57
4106 if(gdrop.and.iep(m).ne.6)ityptl(nptl)=54
4107 if(ghadr)ityptl(nptl)=54
4112 !write(6,'(a,i9)')' ',idptl(nptl)
4119 write(ifch,'(a,5e11.3,$)')' kink:',(pptl(k,nptl),k=1,5)
4120 write(ifch,*)' id: ',idptl(nptl)
4123 c............................no string = resonance...................
4127 if(gdrop)anreso1=anreso1+1
4130 if(nptl.gt.mxptl)call utstop('ProRef: mxptl too small&')
4131 if(iept.eq.0)call idmass(id,am)
4133 pptl(1,nptl)=sngl(ept(1))
4134 pptl(2,nptl)=sngl(ept(2))
4135 pptl(3,nptl)=sngl(ept(3))
4136 pptl(4,nptl)=sngl(ept(4))
4140 if(.not.gdrop)istptl(mm)=41
4142 if(nmes.eq.0.and.nbar.eq.0.and..not.gdrop)ifrptl(1,mm)=nptl
4146 xorptl(1,nptl)=xorptl(1,mm)
4147 xorptl(2,nptl)=xorptl(2,mm)
4148 xorptl(3,nptl)=xorptl(3,mm)
4149 xorptl(4,nptl)=xorptl(4,mm)
4150 tivptl(1,nptl)=xorptl(4,nptl)
4151 call idtau(idptl(nptl),pptl(4,nptl),pptl(5,nptl),taugm)
4152 tivptl(2,nptl)=tivptl(1,nptl)+taugm*(-alog(rangen()))
4155 if(gdrop)ityptl(nptl)=46
4156 if(ghadr)ityptl(nptl)=46
4157 if(iept.eq.6)ityptl(nptl)=47
4158 if(iept.eq.2)ityptl(nptl)=48
4161 if(gdrop)ityptl(nptl)=56
4162 if(ghadr)ityptl(nptl)=56
4163 if(iept.eq.6)ityptl(nptl)=57
4164 if(iept.eq.2)ityptl(nptl)=58
4170 if(ish.ge.3)write(ifch,'(a,5e10.3,i7)')' nucl:'
4171 * ,(pptl(i,nptl),i=1,5),idptl(nptl)
4174 c.......................................................................
4175 c print *,iep(1),iet(1),ityptl(nptl)
4176 1000 call utprix('ProReF',ish,ishini,3)
4177 ctp060829 if(ityptl(nptl).gt.60)print*,ityptl(nptl)
4182 c---------------------------------------------------------------------------------------
4183 subroutine RemoveHadrons(gproj,gtarg,ghadr,m,mm,jcf,icf,ept)
4184 c---------------------------------------------------------------------------------------
4186 include 'epos.incems'
4187 integer jcf(nflav,2),icf(2)
4188 double precision aa(5),ept(5)
4189 logical ghadr,gproj,gtarg
4190 common /cncl/xproj(mamx),yproj(mamx),zproj(mamx)
4191 * ,xtarg(mamx),ytarg(mamx),ztarg(mamx)
4198 call utstop('RemoveHadron : neither proj or targ !&')
4201 call idquacjc(jcf,nqu,naq)
4205 elseif(nqu.gt.naq)then
4207 nbar=(nqu-naq-3)/3 !nbar baryons
4210 nbar=(naq-nqu-3)/3 !nbar antibaryons
4212 if(nmes+nbar.gt.0)ghadr=.true.
4216 !write(ifch,*)'remove meson',mes,' / ',nmes
4217 call gethadron(1,idd,aa,jcf,ept,ir,iret)
4218 call idenco(jcf,icf,iret2)
4219 if(iret.eq.0.and.iret2.eq.0)then
4222 & call utstop('RemoveHadrons: mxptl too small&')
4225 pptl(i,nptl)=sngl(aa(i))
4240 xorptl(1,nptl)=xproj(m)
4241 xorptl(2,nptl)=yproj(m)
4242 xorptl(3,nptl)=zproj(m)
4245 xorptl(1,nptl)=xtarg(m)
4246 xorptl(2,nptl)=ytarg(m)
4247 xorptl(3,nptl)=ztarg(m)
4249 xorptl(4,nptl)=xorptl(4,mm)
4250 tivptl(1,nptl)=xorptl(4,nptl)
4251 call idtau(idptl(nptl),pptl(4,nptl),pptl(5,nptl),taugm)
4252 tivptl(2,nptl)=tivptl(1,nptl)+taugm*(-alog(rangen()))
4255 c deleted: after abstracting a meson,
4256 c check if the NEW remnant is a H-Dibaryon
4259 c remove (anti)baryons
4260 call idquacjc(jcf,nqu,naq)
4263 !write(ifch,*)'remove baryon',nb,' / ',nbar
4265 call gethadron(2,idd,aa,jcf,ept,ir,iret)
4267 call gethadron(3,idd,aa,jcf,ept,ir,iret)
4269 call idenco(jcf,icf,iret2)
4270 if(iret.eq.0.and.iret2.eq.0)then
4273 & call utstop('RemoveHadron: mxptl too small&')
4276 pptl(i,nptl)=sngl(aa(i))
4280 if(nmes.eq.0.and.nb.eq.1)then
4291 xorptl(1,nptl)=xproj(m)
4292 xorptl(2,nptl)=yproj(m)
4293 xorptl(3,nptl)=zproj(m)
4296 xorptl(1,nptl)=xtarg(m)
4297 xorptl(2,nptl)=ytarg(m)
4298 xorptl(3,nptl)=ztarg(m)
4300 xorptl(4,nptl)=xorptl(4,mm)
4301 tivptl(1,nptl)=xorptl(4,nptl)
4302 call idtau(idptl(nptl),pptl(4,nptl),pptl(5,nptl),taugm)
4303 tivptl(2,nptl)=tivptl(1,nptl)+taugm*(-alog(rangen()))
4306 c deleted: after abstracting a (anti)baryon,
4307 c check if the NEW remnant is a H-Dibaryon
4312 c------------------------------------------------------------------
4313 subroutine gethadron(imb,idf,a,jc,ep,ir,iret)
4314 c------------------------------------------------------------------
4315 c goal: emit a hadron (imb= 1 meson, 2 baryon, 3 antibaryon)
4316 c update the remnant flavour and 5-momentum
4318 c idf ,a : hadron id and 5-momentum
4319 c ir : 1 projectile, -1 target remnant
4320 c jc, ep : remnant flavor and 5-momentum
4321 c iret : in case of error, keep correct momentum in remnant
4322 c and lose the quarks of the (not) emitted hadron
4323 c-----------------------------------------------------------------
4326 include 'epos.incems'
4328 double precision s,plc
4329 double precision ep(5),a(5),re(5),p1(5)
4330 integer jc(nflav,2),ifh(3)!,ic(2)
4331 common /ems12/iodiba,bidiba ! defaut iodiba=0. if iodiba=1, study H-Dibaryon
4332 double precision ptm,qcm,u(3),utpcmd,ptt,phi,sxini,strmas
4333 & ,ampt2dro,ampt2str,p5sq,amasex,drangen
4335 call utpri('gethad',ish,ishini,5)
4344 write(ifch,*)'remnant flavour and 5-momentum:',jc, ep, ir
4346 !write(*,'(/a,5f8.3)')'p before: ',ep
4354 c get the id and mass of hadron, the remnant jc is updated
4356 if(imb.eq.1)then ! a meson
4357 ifq=idrafl(iclpt,jc,1,'v',iret)
4358 ifa=idrafl(iclpt,jc,2,'v',iret)
4362 idf=-(ifq*10+ifa*100)
4364 call idmass(idf,amss)
4366 elseif(imb.eq.2)then ! a baryon
4368 ifh(ik)=idrafl(iclpt,jc,1,'v',iret)
4370 call neworder(ifh(1),ifh(2),ifh(3))
4371 idf=ifh(1)*1000+ifh(2)*100+ifh(3)*10
4372 if(ifh(1).ne.ifh(2).and.ifh(2).ne.ifh(3)
4373 $ .and.ifh(1).ne.ifh(3)) idf=2130
4374 if(ifh(1).eq.ifh(2).and.ifh(2).eq.ifh(3))idf=idf+1
4375 call idmass(idf,amss)
4377 elseif(imb.eq.3)then ! an antibaryon
4379 ifh(ik)=idrafl(iclpt,jc,2,'v',iret)
4381 call neworder(ifh(1),ifh(2),ifh(3))
4382 idf=ifh(1)*1000+ifh(2)*100+ifh(3)*10
4383 if(ifh(1).ne.ifh(2).and.ifh(2).ne.ifh(3)
4384 $ .and.ifh(1).ne.ifh(3)) idf=2130
4385 if(ifh(1).eq.ifh(2).and.ifh(2).eq.ifh(3))idf=idf+1
4387 call idmass(idf,amss)
4390 if(iret.ne.0)call utstop('Not enough quark in gethad ???&')
4392 c boost remnant in rest frame
4393 if(ish.ge.6) write (ifch,*) 'on-shell check'
4397 p1(5)=(p1(4)-p1(3))*(p1(4)+p1(3))-p1(2)**2-p1(1)**2
4398 if(p1(5).gt.0d0.and.abs(p1(5)-ep(5)*ep(5)).lt.ep(5))then
4401 if(ish.ge.1)write(ifch,*)'Precision problem in gethad, p:',
4402 & (p1(k),k=1,5),ep(5)*ep(5)
4404 p1(4)=sqrt(p1(3)*p1(3)+p1(2)*p1(2)+p1(1)*p1(1)+p1(5)*p1(5))
4411 strmas=dble(2.*utamnz(jc,4))
4419 !write(ifch,*)'nredo.gt.20 -> only drop'
4420 if(ish.ge.4)write(ifch,*)
4421 & 'Pb with hadron momentum in Gethad !'
4428 ptt=dble(ranpt()*alpdro(2))**2 !pt+pl
4429 if(ptt.ge.sxini)goto 777
4430 sxini=sqrt(sxini-ptt)
4433 ampt2dro=amasex**2d0
4434 ampt2str=strmas**2d0
4438 if(re(5).le.strmas)then
4439 if(ish.ge.6)write(ifch,*)
4440 & 'Pb with initial mass in Gethad, retry',ir
4441 & ,amasex,strmas,sxini,ptm,ptt
4447 if(ish.ge.6)write(ifch,*)'2 body decay',ptm,a(5),re(5)
4448 qcm=utpcmd(ptm,a(5),re(5))
4449 u(3)=2.d0*drangen(qcm)-1.d0
4450 phi=2.d0*dble(pi)*drangen(u(3))
4451 u(1)=sqrt(1.d0-u(3)**2)*cos(phi)
4452 u(2)=sqrt(1.d0-u(3)**2)*sin(phi)
4453 if(u(3).ge.0d0)then !send always hadron backward
4465 re(4)=sqrt(qcm**2+re(5)**2)
4466 a(4)=sqrt(qcm**2+a(5)**2)
4468 if(ish.ge.6)write(ifch,*)'momentum in rest frame : ',re,a
4473 c boost string in collision frame
4474 call utlob2(-1,p1(1),p1(2),p1(3),p1(4),p1(5)
4475 $ ,re(1),re(2),re(3),re(4),81)
4477 p5sq=(re(4)+re(3))*(re(4)-re(3))-(re(1)*re(1)+re(2)*re(2))
4478 if(p5sq.gt.ampt2str)then
4482 write(ifch,*)'Pb with remnant mass -> retry'
4483 write(ifch,*)' m^2:',p5sq,' m_min^2:',ampt2str
4484 write(ifch,*)' momentum four vector:',(re(ii),ii=1,4)
4491 c boost droplet in collision frame
4492 call utlob2(-1,p1(1),p1(2),p1(3),p1(4),p1(5)
4493 $ ,a(1),a(2),a(3),a(4),82)
4495 p5sq=(a(4)+a(3))*(a(4)-a(3))-(a(1)**2.d0+a(2)**2.d0)
4496 if(abs(p5sq-ampt2dro).le.0.1)then
4500 write(ifch,*)'Pb with hadron mass -> retry'
4501 write(ifch,*)' m^2:',p5sq,' m_min^2:',ampt2dro
4502 write(ifch,*)' momentum four vector:',(a(ii),ii=1,4)
4511 if(iret.eq.1)then !If problem with momenta do not update remnant
4514 * write(ifch,*)'no hadron emission in gethad'
4516 else !update the 3-momentum and energy of remnant: ep
4518 if(ish.ge.1.and.abs(ep(4)-re(4)-a(4)).gt.1.e-2*ep(4))then
4519 write(ifmt,*)'Pb with energy conservation in gethad'
4521 write(ifch,*)'Pb with energy conservation :'
4522 write(ifch,*)' p1_ini:',ep(1),' p1:',re(1)+a(1)
4523 write(ifch,*)' p2_ini:',ep(2),' p2:',re(2)+a(2)
4524 write(ifch,*)' p3_ini:',ep(3),' p3:',re(3)+a(3)
4532 write(ifch,*)'get hadron with id and 5-momentum:',idf, a
4539 !write(*,'(a,5f8.3,i5)')'p after: ',sm,iret
4542 c if(abs((a(4)+a(3))*(a(4)-a(3))
4543 c $ -a(2)**2-a(1)**2-a(5)**2).gt.0.3
4544 c $ .and. abs(1.-abs(a(3))/a(4)).gt.0.01)print*,iret,dd
4546 c$$$ if(iodiba.eq.1)then ! for H-dibaryon study ??????????
4547 c$$$ call idenco(jc,ic,iret)
4548 c$$$ if(ic(1).eq.222000.and.ic(2).eq.0)ep(5)=ep(5)-bidiba
4552 write(ifch,*)'new remnant flavour and 5-momentum:',jc, ep
4554 c write(ifmt,*)'get hadron with id and 5-momentum:',idf, a
4555 c write(ifmt,*)'new remnant flavour and 5-momentum:',jc, ep
4557 call utprix('gethad',ish,ishini,5)
4564 c------------------------------------------------------------------
4565 subroutine getdroplet(ir,ic,jc,ep,a,pass,mdrmax)
4566 c------------------------------------------------------------------
4567 c emit a droplet, update the remnant string flavour and 5-momentum
4570 c ir ........ 1 projectile, -1 target remnant
4571 c ep ........ remnant 5-momentum
4572 c jc ........ remnant jc
4574 c pass ... .true. = successful droplet emission
4575 c ic, ep ....... droplet ic and 5-momentum
4576 c jc, a ........ remnant string jc and 5-momentum
4577 c .false. = unsuccessful
4578 c jc, ep .... unchanged,
4579 c considered as droplet jc and 5-momentum
4580 c-----------------------------------------------------------------
4583 include 'epos.incems'
4584 double precision ep(5),a(5),p1(5),re(5),eps,amasex,mdrmax
4585 double precision xxx,rr,alp,p5sq,xmin,xmax,ampt2str
4586 & ,sxini,strmas,xxxmax,xxxmin,ampt2dro,mdrmaxi
4587 parameter(eps=1.d-20)
4588 integer jc(nflav,2),jcini(nflav,2),jcfin(nflav,2),ifh(3),ic(2)
4589 integer icx(2) !,icxx(2)
4592 double precision s,plc,ptm,qcm,u(3),utpcmd,ptt,drangen,phi
4594 call utpri('getdro',ish,ishini,4)
4609 call idquacjc(jc,nqu,naq)
4624 write(ifch,*)'remnant flavour and 5-momentum:',jc,ep,npart
4627 c get id of string ends, the remnant string jc is updated
4629 if(npart.lt.3.and.ep(5).lt.mdrmax.and.iclpt.ne.4)then !light droplet with few quarks
4632 elseif(npart.lt.3)then !few quarks but heavy, add some quarks to extract a q-qbar string (should not exit directly because of the large mass)
4633 ifq=jdrafl(iclpt,jcini,2,iret2)
4634 if(nqu.eq.1.and.naq.eq.1)then
4640 call utstop('This should not happen (getdrop) !&')
4642 elseif((nqu.eq.2.and.naq.le.2).or.(nqu.le.2.and.naq.eq.2))then
4645 elseif(naq.eq.0)then
4648 elseif(nqu.eq.0)then
4651 else !There is enough q or aq to do qq-q string
4654 if(jcini(4,1)-jcini(4,2).eq.0)then !if c-cbar
4661 c One chooses the first q or aq
4665 if(jcini(4,1)+jcini(4,2).ne.0)then !if some charm take it out
4666 if(jcini(4,1).ne.0)then
4673 elseif(rrr.gt.float(naq)/float(npart))then
4681 c One chooses the second one
4685 if(idps.eq.1.and.jcini(4,1).ne.0)then !if some charm take it out
4687 elseif(idms.eq.1.and.jcini(4,2).ne.0)then !if some charm take it out
4689 elseif(rrr.gt.float(naq)/float(npart))then
4690 if(idps.eq.1.and.nqu.ge.2)then
4696 if(idms.eq.1.and.naq.ge.2)then
4703 c If there is already 2 q or 2 aq as string end, we know that we need
4704 c a third one to complete the string
4708 if(idps.eq.1.and.idms.ne.5)idms=1
4709 if(idms.eq.1.and.idps.ne.5)idps=1
4716 write(ifch,*)'remnant string ends :',idps,idms
4719 if(idps.ne.5.and.idms.ne.5)then ! q-aq string
4720 if(jcini(4,1).eq.1)then
4721 ifq=idrafl(iclpt,jcini,1,'c',iret)
4723 ifq=idrafl(iclpt,jcini,1,'v',iret)
4725 if(jcini(4,1).eq.1)then
4726 ifa=idrafl(iclpt,jcini,2,'c',iret)
4728 ifa=idrafl(iclpt,jcini,2,'v',iret)
4733 elseif(idps.eq.5)then ! qq-q string
4735 if(jcini(4,1).ne.0)then
4736 ifh(ik)=idrafl(iclpt,jcini,1,'c',iret)
4738 ifh(ik)=idrafl(iclpt,jcini,1,'v',iret)
4740 jcfin(ifh(ik),1)=jcfin(ifh(ik),1)+1
4743 elseif(idms.eq.5)then !aqaq-aq string
4745 if(jcini(4,2).ne.0)then
4746 ifh(ik)=idrafl(iclpt,jcini,2,'c',iret)
4748 ifh(ik)=idrafl(iclpt,jcini,2,'v',iret)
4750 jcfin(ifh(ik),2)=jcfin(ifh(ik),2)+1
4754 if(iret.ne.0)call utstop('Not enough quark in getdro ???&')
4755 if(jcini(4,1)+jcini(4,2).ne.0)
4756 & call utstop('There is sitll charm quark in getdro???&')
4760 call idenco(jcini,icx,iret)
4762 call utstop('Exotic flavor in getdroplet !&')
4765 idx=idtra(icx,0,0,3)
4766 if(idx.ne.0)call idmass(idx,amx)
4770 c boost remnant in rest frame
4771 if(ish.ge.6) write (ifch,*) 'on-shell check'
4775 p1(5)=(p1(4)-p1(3))*(p1(4)+p1(3))-p1(2)**2-p1(1)**2
4776 if(p1(5).gt.0d0.and.abs(p1(5)-ep(5)*ep(5)).lt.ep(5))then
4779 if(ish.ge.2)write(ifch,*)'Precision problem in getdro, p:',
4780 & (p1(k),k=1,5),ep(5)*ep(5)
4782 p1(4)=sqrt(p1(3)*p1(3)+p1(2)*p1(2)+p1(1)*p1(1)+p1(5)*p1(5))
4784 if(ish.ge.6) write (ifch,*) 'boost vector:',p1
4786 c limits for momenta
4793 amasex=dble(fad*utamnz(jcini,mamod))
4794 strmas=dble(fas*utamnz(jcfin,mamos))
4803 amasex=dble(utamnz(jcini,mamod))
4804 strmas=dble(utamnz(jcfin,mamos))
4805 elseif(nredo.gt.20)then
4806 !write(ifch,*)'nredo.gt.20 -> only drop'
4807 if(ish.ge.4)write(ifch,*)
4808 & 'Pb with string mass in Getdrop, continue with gethad'
4816 ptt=dble(ranpt()*alpdro(2))**2 !pt+pl
4817 if(ptt.ge.sxini)goto 777
4818 sxini=sqrt(sxini-ptt)
4821 ampt2dro=amasex**2d0
4822 ampt2str=strmas**2d0
4823 if(ampt2dro.gt.mdrmaxi)then
4824 mdrmaxi=2d0*ampt2dro
4825 c write(ifmt,*)'Warning Mmin>Mmax in Getdroplet'
4828 xxxmax=min(mdrmaxi,(sxini-strmas)**2) !strmas/(strmas+ampt2)
4831 if(xxxmin.gt.xxxmax)then
4832 !write(ifch,*)'Warning Mmin>sxini -> only drop'
4833 if(ish.ge.4)write(ifch,*)
4834 & 'Pb with ampt2 in Getdrop, retry',nredo,ir
4835 & ,ampt2dro,ampt2str,xxxmin,xxxmax,sxini,ptt,mdrmaxi
4847 if(dabs(alp-1.d0).lt.eps)then
4848 xxx=xmax**rr*xmin**(1d0-rr)
4850 xxx=(rr*xmax**(1d0-alp)+(1d0-rr)*xmin**(1d0-alp))
4854 c write(ifch,*)'ini',xmin,xxx,xmax,rr,ampt2dro
4855 c & ,(sxini-sqrt(xxx)),ampt2str,p1(5)
4861 if(a(5).le.strmas)then
4862 if(ish.ge.6)write(ifch,*)
4863 & 'Pb with initial mass in Getdrop, retry',ir
4864 & ,xmin,xxx,xmax,rr,ampt2dro,ampt2str
4870 if(ish.ge.6)write(ifch,*)'2 body decay',ptm,re(5),a(5)
4871 qcm=utpcmd(ptm,re(5),a(5))
4872 u(3)=2.d0*drangen(qcm)-1.d0
4873 phi=2.d0*dble(pi)*drangen(u(3))
4874 u(1)=sqrt(1.d0-u(3)**2)*cos(phi)
4875 u(2)=sqrt(1.d0-u(3)**2)*sin(phi)
4876 if(u(3).lt.0d0)then !send always droplet backward
4888 re(4)=sqrt(qcm**2+re(5)**2)
4889 a(4)=sqrt(qcm**2+a(5)**2)
4891 if(ish.ge.6)write(ifch,*)'momentum in rest frame : ',re,a
4897 c boost string in collision frame
4898 call utlob2(-1,p1(1),p1(2),p1(3),p1(4),p1(5)
4899 $ ,a(1),a(2),a(3),a(4),71)
4901 p5sq=(a(4)+a(3))*(a(4)-a(3))-(a(1)**2.d0+a(2)**2.d0)
4902 if(p5sq.gt.ampt2str)then
4906 write(ifch,*)'Pb with string mass -> retry'
4907 write(ifch,*)' m^2:',p5sq,' m_min^2:',ampt2str
4908 write(ifch,*)' momentum four vector:',(a(ii),ii=1,4)
4915 c boost droplet in collision frame
4916 call utlob2(-1,p1(1),p1(2),p1(3),p1(4),p1(5)
4917 $ ,re(1),re(2),re(3),re(4),72)
4919 p5sq=(re(4)+re(3))*(re(4)-re(3))-(re(1)*re(1)+re(2)*re(2))
4920 if(p5sq.gt.ampt2dro)then
4924 write(ifch,*)'Pb with droplet mass -> retry'
4925 write(ifch,*)' m^2:',p5sq,' m_min^2:',ampt2dro
4926 write(ifch,*)' momentum four vector:',(re(ii),ii=1,4)
4932 if(ish.ge.1.and.abs(ep(4)-re(4)-a(4)).gt.1.e-2*ep(4))then
4933 write(ifmt,*)'Pb with energy conservation in getdro'
4935 write(ifch,*)'Pb with energy conservation :'
4936 write(ifch,*)' p1_ini:',ep(1),' p1:',re(1)+a(1)
4937 write(ifch,*)' p2_ini:',ep(2),' p2:',re(2)+a(2)
4938 write(ifch,*)' p3_ini:',ep(3),' p3:',re(3)+a(3)
4942 c If OK, save flavours of droplet and string
4954 write(ifch,*)'droplet:'
4957 write(ifch,*)'string remnant:'
4963 call utprix('getdro',ish,ishini,4)
4966 c-----------------------------------------------------
4967 subroutine neworder(n1, n2, n3)
4968 c-----------------------------------------------------
4969 c make 3 integers ordered like 1 2 3
4970 c------------------------------------------------------
4981 elseif(n3.lt.n2)then
4988 c-----------------------------------------------------------------------
4990 c-----------------------------------------------------------------------
4991 c transforms ic to id such that only hadrons have nonzero id
4992 c-----------------------------------------------------------------------
4994 integer idt(3,nidt),ic(2)
4996 * 100000,100000, 110 ,100000,010000, 120 ,010000,010000, 220
4997 *,100000,001000, 130 ,010000,001000, 230 ,001000,001000, 330
4998 *,100000,000100, 140 ,010000,000100, 240 ,001000,000100, 340
4999 *,000100,000100, 440
5000 *,300000,000000,1111 ,210000,000000,1120 ,120000,000000,1220
5001 *,030000,000000,2221 ,201000,000000,1130 ,111000,000000,1230
5002 *,021000,000000,2230 ,102000,000000,1330 ,012000,000000,2330
5003 *,003000,000000,3331 ,200100,000000,1140 ,110100,000000,1240
5004 *,020100,000000,2240 ,101100,000000,1340 ,011100,000000,2340
5005 *,002100,000000,3340 ,100200,000000,1440 ,010200,000000,2440
5006 *,001200,000000,3440 ,000300,000000,4441/
5009 if(ic(1).eq.0.and.ic(2).eq.0)then
5010 if(rangen().ge.0.5)then
5022 if(ic(2).eq.idt(1,i).and.ic(1).eq.idt(2,i))idtr2=-idt(3,i)
5023 if(ic(1).eq.idt(1,i).and.ic(2).eq.idt(2,i))idtr2=idt(3,i)
5028 c----------------------------------------------------------------------
5029 subroutine emsini(e,idpj,idtg)
5030 c----------------------------------------------------------------------
5031 c energy-momentum sharing initializations
5032 c----------------------------------------------------------------------
5034 include 'epos.incems'
5035 include 'epos.incsem'
5036 common/cemsr5/at(0:1,0:6)
5038 common/cems10/a(0:ntypmx),b(0:ntypmx),d(0:ntypmx)
5039 common/ems6/ivp0,iap0,idp0,isp0,ivt0,iat0,idt0,ist0
5040 double precision d,a,b,plc,s,amd,dcel,xvpr,xdm,at,xdm2
5042 common/cems13/xvpr(0:3)
5051 c alpha (0=0, 1=s, 2=v, 4=d, 8=f)
5064 c beta (0=0, 1=s, 2=v, 4=d, 8=f)
5078 c alpha_trailing and beta_trailing (0=meson, 1=baryon;
5079 c 0=no excit, 1=nondiffr, 2=diffr,
5080 c 6=nondiffr but no pomeron)
5083 at(0,1)=dble(alpndi)
5085 at(0,3)=dble(alpndi)
5086 at(0,6)=dble(alpndi)
5088 at(1,1)=dble(alpndi)
5090 at(1,3)=dble(alpndi)
5091 at(1,6)=dble(alpndi)
5093 c minimal string masses ( i+j, each one: 0=0, 1=s, 2=v, 4=d, 8=f)
5097 ammn(2)=dble(ammsqq)+amd
5098 ammn(3)=dble(ammsqq)+amd
5099 ammn(4)=dble(ammsqq)+amd
5100 ammn(5)=dble(ammsqd)+amd
5101 ammn(6)=dble(ammsqd)+amd
5103 ammn(8)=dble(ammsdd)+amd
5104 ammn(9)=dble(ammsqd)+amd
5105 ammn(10)=dble(ammsqd)+amd
5106 ammn(12)=dble(ammsqd)+amd
5109 c minimal pomeron masses (0=0, 1=softPom, 2=regge, 3=hard)
5114 amprmn(3)=dsqrt(4d0*dble(q2min))
5116 c cutoff for virtual pomeron (0=0, 1=soft Pom, 2=regge, 3=hard)
5119 xvpr(1)=dble(cumpom**2)/s
5120 xvpr(2)=dble(cumpom**2)/s
5123 c minimal remnant masses (0=meson, 1=baryon)
5126 call idmass(idpj,ampj)
5127 if(iabs(idpj).gt.1000)then
5129 ampmn(1)=dble(ampj)+xdm
5131 ampmn(0)=dble(ampj)+xdm
5134 call idmass(idtg,amtg)
5135 if(iabs(idtg).gt.1000)then
5137 amtmn(1)=dble(amtg)+xdm
5139 amtmn(0)=dble(amtg)+xdm
5143 c minimal excitation masses (0=meson, 1=baryon
5144 c 0=no excit, 1=nondiffr, 2=diffr,
5145 c 6=nondiffr but no pomeron)
5148 c to take into account increase of mean pt in inelastic remnants
5149 c if(isplit.eq.1)xdm=max(2d0*max(0d0,sqrt(log(s))-2.5d0),xdm)
5151 amemn(0,1)=xdm2+0.31d0
5153 amemn(0,2)=xdm2+0.31d0
5155 amemn(0,3)=xdm2+0.31d0
5156 c amemn(0,6)=xdm2+0.31d0
5159 amemn(1,1)=xdm2+0.15d0 !+0.15d0
5161 amemn(1,2)=xdm2+0.15d0
5163 amemn(1,3)=xdm2+0.15d0
5164 c amemn(1,6)=xdm2+0.15d0
5167 c maximal excitation masses (0=no excit, 1=nondiffr, 2=diffr)
5173 if(idpj.gt.1000)then ! baryon
5175 c initial quark configuration
5181 elseif(idpj.lt.-1000)then ! antibaryon
5183 c initial quark configuration
5191 c initial quark configuration
5199 if(idtg.gt.1000)then ! baryon
5201 c initial quark configuration
5207 elseif(idtg.lt.-1000)then ! antibaryon
5209 c initial quark configuration
5217 c initial quark configuration
5226 c eikonal parameters
5228 dcel=dble(chad(iclpro)*chad(icltar))
5254 c-----------------------------------------------------------------------
5256 c-----------------------------------------------------------------------
5258 c-----------------------------------------------------------------------
5261 include 'epos.incems'
5263 common/emsptl/nppr(npommx,kollmx),npproj(mamx),nptarg(mamx)
5265 call utpri('emsigr',ish,ishini,5)
5267 do k=1,koll !----k-loop---->
5269 c determine length of k-th line of grid
5271 o=max(1.e-5,min(sngl(om1intc(k)),float(npommx)))!if GFF used for propo
5272 if(ish.ge.7)write(ifch,*)'emsigr:k,o',k,o
5281 if(ish.ge.7)write(ifch,*)'emsigr:n,p',n,p
5282 if((p.gt.1e-4.or.n.lt.int(o)).and.n.lt.npommx
5283 *.and.n.lt.nprmax)goto 10
5285 if(ish.ge.5)write(ifch,*)'emsigr:nmax,b',n,bk(k)
5295 c initial value for interaction type
5340 enddo ! <----k-loop-----
5342 call utprix('emsigr',ish,ishini,5)
5346 c-----------------------------------------------------------------------
5348 c-----------------------------------------------------------------------
5349 c initialize projectile and target
5350 c-----------------------------------------------------------------------
5353 include 'epos.incems'
5355 common/emsptl/nppr(npommx,kollmx),npproj(mamx),nptarg(mamx)
5358 common/ems6/ivp0,iap0,idp0,isp0,ivt0,iat0,idt0,ist0
5359 common /cncl/xproj(mamx),yproj(mamx),zproj(mamx)
5360 * ,xtarg(mamx),ytarg(mamx),ztarg(mamx)
5362 double precision dcel,s,plc
5364 c initialize projectile
5378 xpmn(i)=(amemn(idp(i),0)+ampmn(idp(i)))**2/s
5379 xpmx(i)=dmin1(1d0,(amemx(0)+ampmn(idp(i)))**2/s)
5380 xpos(i)=0.9d0*(amemx(0)+ampmn(idp(i)))**2/s
5381 xppmx(i)=0.5d0/(1d0+1d0/dble(maproj)**0.3d0)!1d0-dsqrt(xpmn(i))/maproj
5382 xmpmx(i)=0.5d0/(1d0+1d0/dble(matarg)**0.3d0)!1d0-dsqrt(xpmn(i))/matarg
5383 xmpmn(i)=xpmn(i)/xppmx(i)
5384 xppmn(i)=xpmn(i)/xmpmx(i)
5406 xtmn(j)=(amemn(idt(j),0)+amtmn(idt(j)))**2/s
5407 xtmx(j)=dmin1(1d0,(amemx(0)+amtmn(idt(j)))**2/s)
5408 xtos(j)=0.9d0*(amemx(0)+amtmn(idt(j)))**2/s
5409 xmtmx(j)=0.5d0/(1d0+1d0/dble(matarg)**0.3d0)!1d0-dsqrt(xtmn(j))/matarg
5410 xptmx(j)=0.5d0/(1d0+1d0/dble(maproj)**0.3d0)!1d0-dsqrt(xtmn(j))/maproj
5411 xptmn(j)=xtmn(j)/xmtmx(j)
5412 xmtmn(j)=xtmn(j)/xptmx(j)
5424 c-----------------------------------------------------------------------
5426 c-----------------------------------------------------------------------
5427 c completes /cptl/ for nucleons, checks for no interaction
5429 c-----------------------------------------------------------------------
5431 include 'epos.incems'
5432 common/nucl3/phi,bimp
5433 common/col3/ncol,kolpt
5434 integer kolpz(mamx),koltz(mamx)
5436 call utpri('emszz ',ish,ishini,6)
5441 if(iokoll.eq.1)then ! precisely matarg collisions
5453 if(ish.ge.7)write(ifch,*)'k,itpr,ncol,ncolx',k,itpr(k),ncol,ncolx
5454 if(itpr(k).eq.0)goto 8
5455 if(itpr(k).eq.1.or.itpr(k).eq.3)ncoli=ncoli+1
5461 tivptl(2,i)=coord(4,k)
5464 tivptl(2,maproj+j)=coord(4,k)
5466 if(ncolx.ne.ncol)write(6,*)'ncolx,ncol:', ncolx,ncol
5467 if(ncolx.ne.ncol)call utstop('********ncolx.ne.ncol********&')
5468 if(ncol.eq.0)goto1001
5470 c determine npj, ntg
5479 if(itpr(k).gt.0)then
5482 kolpz(ip)=kolpz(ip)+1
5483 koltz(it)=koltz(it)+1
5488 if(kolpz(ip).gt.0)npj=npj+1
5492 if(koltz(it).gt.0)ntg=ntg+1
5494 c write(6,*)'npj,ntg,npj+ntg:',npj,ntg,npj+ntg
5510 !print*,' ===== ',kolevt,koievt,' ====='
5517 write(ifch,115)iorptl(n),jorptl(n),n,istptl(n)
5518 *,tivptl(1,n),tivptl(2,n)
5520 115 format(1x,'/cptl/',2i6,2i10,2(e10.3,1x))
5524 call utprix('emszz ',ish,ishini,6)
5530 write(ifch,*)' ***** no interaction!!!'
5531 write(ifch,*)' ***** ncol=0 detected in emszz'
5538 c-----------------------------------------------------------------------
5539 subroutine ProCop(i,ii)
5540 c-----------------------------------------------------------------------
5541 c Propose Coordinates of remnants from active projectile nucleons
5542 c-----------------------------------------------------------------------
5544 include 'epos.incems'
5547 double precision xmptmp,aproj
5549 common/emsptl/nppr(npommx,kollmx),npproj(mamx),nptarg(mamx)
5551 double precision s,plc
5555 idptl(nptl)=idptl(ii)*100+99 !100*10**idp(i)+iep(i)
5567 if(lproj(i).gt.1)then
5573 if(itpr(k).gt.0.and.z.gt.zmax)then
5581 c if(kolz.eq.0)call utstop(' kolz=0 (proj)&')
5588 xorptl(1,nptl)=xorptl(1,ii)
5589 xorptl(2,nptl)=xorptl(2,ii)
5590 xorptl(3,nptl)=xorptl(3,ii)
5595 icrmn(1)=icproj(1,i)
5596 icrmn(2)=icproj(2,i)
5597 aproj=dble(max(amproj,fremnux2(icrmn)))
5598 c aprojex=max(ampmn(idp(i))+amemn(idp(i),iep(i))
5599 c & ,dble(fremnux(icrmn)))
5600 xmptmp=(aproj**2+xxp(i)*xxp(i)+xyp(i)*xyp(i))
5602 if(iep(i).eq.6)xmptmp=max(xmptmp,xmp(i))
5603 xpos(i)=1.1d0*xpp(i)*xmptmp
5604 if(xmptmp.gt.1.d0)then
5606 if(ish.ge.1)write(ifmt,*)'Warning in ProCop, Remnant mass too low'
5609 pptl(1,nptl)=sngl(xxp(i))
5610 pptl(2,nptl)=sngl(xyp(i))
5611 pptl(3,nptl)=sngl((xpp(i)-xmptmp)*plc/2d0)
5612 pptl(4,nptl)=sngl((xpp(i)+xmptmp)*plc/2d0)
5615 c write(ifmt,*)'ProCop',i,nptl
5621 c-----------------------------------------------------------------------
5622 subroutine ProCot(j,jj)
5623 c-----------------------------------------------------------------------
5624 c Propose Coordinates of remnants from active targets nucleons
5625 c-----------------------------------------------------------------------
5628 include 'epos.incems'
5630 double precision xpttmp,atarg
5632 common/emsptl/nppr(npommx,kollmx),npproj(mamx),nptarg(mamx)
5634 double precision s,plc
5639 idptl(nptl)=idptl(jj)*100+99 !100*10**idt(j)+iet(j)
5651 if(ltarg(j).gt.1)then
5657 if(itpr(k).gt.0.and.z.lt.zmin)then
5665 c if(kolz.eq.0)call utstop(' kolz=0 (targ)&')
5672 xorptl(1,nptl)=xorptl(1,jj)
5673 xorptl(2,nptl)=xorptl(2,jj)
5674 xorptl(3,nptl)=xorptl(3,jj)
5679 icrmn(1)=ictarg(1,j)
5680 icrmn(2)=ictarg(2,j)
5681 atarg=dble(max(amtarg,fremnux2(icrmn)))
5682 c atargex=max(amtmn(idt(j))+amemn(idt(j),iet(j))
5683 c & ,dble(fremnux(icrmn)))
5684 xpttmp=(atarg**2+xxt(j)*xxt(j)+xyt(j)*xyt(j))
5686 if(iet(j).eq.6)xpttmp=max(xpttmp,xpt(j))
5687 xtos(j)=1.1d0*xpttmp*xmt(j)
5688 if(xpttmp.gt.1.d0)then
5690 if(ish.ge.1)write(ifch,*)'Warning in ProCot, Remnant mass too low'
5693 pptl(1,nptl)=sngl(xxt(j))
5694 pptl(2,nptl)=sngl(xyt(j))
5695 pptl(3,nptl)=sngl((xpttmp-xmt(j))*plc/2d0)
5696 pptl(4,nptl)=sngl((xpttmp+xmt(j))*plc/2d0)
5699 c write(ifmt,*)'ProCot',j,nptl
5704 c-----------------------------------------------------------------------
5705 subroutine emswrp(i,ii)
5706 c-----------------------------------------------------------------------
5708 include 'epos.incems'
5711 double precision p5sq
5713 common/emsptl/nppr(npommx,kollmx),npproj(mamx),nptarg(mamx)
5714 double precision s,plc
5715 parameter(eps=1.e-5)
5717 if(npproj(i).eq.0)then
5718 write(*,*)'emswrp i ii',i,ii
5719 call utstop('emswrp with npproj=0 should never happen !&')
5721 c t=xorptl(4,kolp(i))
5727 c idptl(nptl)=idptl(ii)*100+99 !100*10**idp(i)+iep(i)
5731 c jorptl(nptl)=kolp(i)
5734 c xorptl(1,nptl)=xorptl(1,ii)
5735 c xorptl(2,nptl)=xorptl(2,ii)
5736 c xorptl(3,nptl)=xorptl(3,ii)
5745 pptl(1,mm)=sngl(xxp(i))
5746 pptl(2,mm)=sngl(xyp(i))
5747 pptl(3,mm)=sngl((xpp(i)-xmp(i))*plc/2d0)
5748 pptl(4,mm)=sngl((xpp(i)+xmp(i))*plc/2d0)
5749 if(pptl(4,mm).lt.-eps)call utstop('E pro<0 !&')
5750 p5sq=xpp(i)*plc*xmp(i)*plc-xxp(i)*xxp(i)-xyp(i)*xyp(i)
5751 if(p5sq.gt.1.d-10)then
5752 pptl(5,mm)=sngl(dsqrt(p5sq))
5755 write(ifch,*)'problem with mass for projectile, '
5756 & ,'continue with zero mass'
5757 write(ifch,*)i,mm,xxp(i),xyp(i),xpp(i),xmp(i),p5sq
5770 c-----------------------------------------------------------------------
5771 subroutine emswrt(j,jj)
5772 c-----------------------------------------------------------------------
5775 include 'epos.incems'
5777 double precision p5sq
5779 common/emsptl/nppr(npommx,kollmx),npproj(mamx),nptarg(mamx)
5780 double precision s,plc
5781 parameter(eps=1.e-5)
5783 if(nptarg(j).eq.0)then
5785 write(*,*)'emswrt j jj',j,jj
5786 call utstop('emswrt with nptarg=0 should never happen !&')
5788 c t=xorptl(4,kolt(j))
5794 c idptl(nptl)=idptl(jj)*100+99 !100*10**idp(i)+iep(i)
5798 c jorptl(nptl)=kolt(j)
5801 c xorptl(1,nptl)=xorptl(1,jj)
5802 c xorptl(2,nptl)=xorptl(2,jj)
5803 c xorptl(3,nptl)=xorptl(3,jj)
5812 pptl(1,mm)=sngl(xxt(j))
5813 pptl(2,mm)=sngl(xyt(j))
5814 pptl(3,mm)=sngl((xpt(j)-xmt(j))*plc/2d0)
5815 pptl(4,mm)=sngl((xpt(j)+xmt(j))*plc/2d0)
5816 if(pptl(4,mm).lt.-eps)call utstop('E targ<0 !&')
5817 p5sq=xpt(j)*plc*xmt(j)*plc-xxt(j)*xxt(j)-xyt(j)*xyt(j)
5818 if(p5sq.gt.1.d-10)then
5819 pptl(5,mm)=sngl(dsqrt(p5sq))
5822 write(ifch,*)'problem with mass for target, '
5823 & ,'continue with zero mass'
5824 write(ifch,*)j,mm,xxt(j),xyt(j),xpt(j),xmt(j),p5sq
5836 c-----------------------------------------------------------------------
5837 subroutine emswrpom(k,i,j)
5838 c-----------------------------------------------------------------------
5841 include 'epos.incems'
5844 common/emsptl/nppr(npommx,kollmx),npproj(mamx),nptarg(mamx)
5845 double precision s,px,py,plc
5848 if(idpr(n,k).eq.0.or.ivpr(n,k).eq.0)goto30
5851 px=xxp1pr(n,k)+xxp2pr(n,k)+xxm1pr(n,k)+xxm2pr(n,k)
5852 py=xyp1pr(n,k)+xyp2pr(n,k)+xym1pr(n,k)+xym2pr(n,k)
5853 pptl(1,nptl)=sngl(px)
5854 pptl(2,nptl)=sngl(py)
5855 pptl(3,nptl)=sngl(dsqrt(xpr(n,k))*dsinh(ypr(n,k))*plc)
5856 pptl(4,nptl)=sngl(dsqrt(xpr(n,k))*dcosh(ypr(n,k))*plc)
5857 pptl(5,nptl)=sngl(dsqrt(xpr(n,k)*plc*plc-px*px-py*py))
5858 ! print*,pptl(5,nptl)/plc
5859 idptl(nptl)=idpr(n,k)*10000
5864 idptl(nptl)=idptl(nptl)*100+99
5870 xorptl(1,nptl)=coord(1,k)
5871 xorptl(2,nptl)=coord(2,k)
5872 xorptl(3,nptl)=coord(3,k)
5873 xorptl(4,nptl)=coord(4,k)
5874 tivptl(1,nptl)=coord(4,k)
5875 tivptl(2,nptl)=coord(4,k)
5876 if(idpr(n,k).eq.1)then
5878 elseif(idpr(n,k).eq.2)then
5880 elseif(idpr(n,k).eq.3)then
5883 call utstop('emswrpom: unknown id&')
5893 cc--------------------------------------------------------------------------
5894 c subroutine reaction(idpj,idtg,ireac)
5895 cc--------------------------------------------------------------------------
5896 cc returns reaction code ireac
5897 cc--------------------------------------------------------------------------
5902 c call idchrg(idpj,cp)
5903 c call idchrg(idtg,ct)
5905 c if(iap.gt.100)then
5906 c if(iat.gt.100)then
5920 c elseif(iat.eq.11.or.iat.eq.12.or.iat.eq.22)then
5933 c elseif(iap.eq.11.or.iap.eq.12.or.iap.eq.22)then
5934 c if(iat.gt.100)then
5940 c elseif(iat.eq.11.or.iat.eq.12.or.iat.eq.22)then
5946 c if(iat.gt.100)then
5952 c elseif(iat.eq.11.or.iat.eq.12.or.iat.eq.22)then
5961 c-----------------------------------------------------------------------
5962 subroutine xEmsI1(iii,kc,omlog)
5963 c-----------------------------------------------------------------------
5964 c plot omlog vs iter
5965 c plot nr of pomerons vs iter
5966 c plot number of collisions vs iter
5967 c-----------------------------------------------------------------------
5970 include 'epos.incems'
5971 include 'epos.incsem'
5974 common/cmc/ot(0:nbin),zz(0:nbin),i(0:nbin)
5975 *,yt1,yt2,kx(0:nbin)
5977 common/cmc1/xp(0:nbim),xt(0:nbim),x(0:nbim),o(0:nbim)
5980 double precision xp,xt,x,omlog,om1intbc
5982 double precision plc,s,seedp
5985 c if(iemsi2.eq.0)call utstop('ERROR in XemsI1: iemsi2 = 0&')
5994 c if(itpr(ko).gt.0)then
5995 if(nprt(ko).gt.0)then
6002 elseif(iii.eq.2)then
6010 ctp060829 ip=iproj(ko)
6011 ctp060829 it=itarg(ko)
6012 om1i=sngl(om1intbc(bk(ko)))
6016 om1g=sngl(om1intbc(bk(ko)))
6019 if(rangen().lt.1.-exp(-om1i))then
6022 if(rangen().lt.1.-exp(-om1g))then
6030 write(ce,'(f8.2)')sngl(plc)
6032 write(ifhi,'(a)') '!##################################'
6033 write(ifhi,'(a,i3)') '! log omega for event ',nrevt+1
6034 write(ifhi,'(a)') '!##################################'
6035 write(ifhi,'(a,i1)') 'openhisto name omega-',nrevt+1
6036 write(ifhi,'(a)') 'htyp lin'
6037 write(ifhi,'(a)') 'xmod lin ymod lin'
6038 write(ifhi,'(a,2e11.3)')'xrange',x1,x2
6039 write(ifhi,'(a)') 'yrange auto auto '
6040 write(ifhi,'(a)') 'text 0 0 "xaxis iteration"'
6041 write(ifhi,'(a)') 'text 0 0 "yaxis ln[W]"'
6042 write(ifhi,'(a,a)') 'text 0.5 0.90 "E ='//ce//'"'
6043 write(ifhi,'(a)') 'array 2'
6045 write(ifhi,'(2e11.3)')float(k),o(k)
6047 write(ifhi,'(a)') ' endarray'
6048 write(ifhi,'(a)') 'closehisto plot 0'
6050 write(ifhi,'(a)') '!##################################'
6051 write(ifhi,'(a,i3)')'! nr of coll`s for event ',nrevt+1
6052 write(ifhi,'(a)') '!##################################'
6053 write(ifhi,'(a,i1)') 'openhisto name coll-',nrevt+1
6054 write(ifhi,'(a)') 'htyp lin'
6055 write(ifhi,'(a)') 'xmod lin ymod lin'
6056 write(ifhi,'(a,2e11.3)')'xrange',x1,x2
6057 write(ifhi,'(a)') 'text 0 0 "xaxis iteration"'
6058 write(ifhi,'(a)') 'text 0 0 "yaxis nr of collisions"'
6059 write(ifhi,'(a)') 'array 2'
6061 write(ifhi,'(2e11.3)')float(k),float(kx(k))
6063 write(ifhi,'(a)') ' endarray'
6064 write(ifhi,'(a)') 'closehisto plot 0-'
6065 write(ifhi,'(a)') 'openhisto'
6066 write(ifhi,'(a)') 'htyp lin'
6067 write(ifhi,'(a)') 'array 2'
6069 write(ifhi,'(2e11.3)')float(k),float(kollx)
6071 write(ifhi,'(a)') ' endarray'
6072 write(ifhi,'(a)') 'closehisto plot 0-'
6073 write(ifhi,'(a)') 'openhisto'
6074 write(ifhi,'(a)') 'htyp lin'
6075 write(ifhi,'(a)') 'array 2'
6077 write(ifhi,'(2e11.3)')float(k),float(kollg)
6079 write(ifhi,'(a)') ' endarray'
6080 write(ifhi,'(a)') 'closehisto plot 0'
6082 write(ifhi,'(a)') '!##################################'
6083 write(ifhi,'(a,i3)')'! nr of pom`s for event ',nrevt+1
6084 write(ifhi,'(a)') '!##################################'
6085 write(ifhi,'(a,i1)') 'openhisto name pom-',nrevt+1
6086 write(ifhi,'(a)') 'htyp lin'
6087 write(ifhi,'(a)') 'xmod lin ymod lin'
6088 write(ifhi,'(a,2e11.3)')'xrange',x1,x2
6089 write(ifhi,'(a)') 'text 0 0 "xaxis iteration"'
6090 write(ifhi,'(a)') 'text 0 0 "yaxis nr of Pomerons"'
6091 write(ifhi,'(a)') 'array 2'
6093 write(ifhi,'(2e11.3)')float(k),zz(k)
6095 write(ifhi,'(a)') ' endarray'
6096 if(sum.lt.4*zz(nbin))then
6097 write(ifhi,'(a)') 'closehisto plot 0-'
6098 write(ifhi,'(a)') 'openhisto'
6099 write(ifhi,'(a)') 'htyp lin'
6100 write(ifhi,'(a)') 'array 2'
6102 write(ifhi,'(2e11.3)')float(k),sum
6104 write(ifhi,'(a)') ' endarray'
6105 write(ifhi,'(a)') 'closehisto plot 0-'
6106 write(ifhi,'(a)') 'openhisto'
6107 write(ifhi,'(a)') 'htyp lin'
6108 write(ifhi,'(a)') 'array 2'
6110 write(ifhi,'(2e11.3)')float(k),sumg
6112 write(ifhi,'(a)') ' endarray'
6114 write(ifhi,'(a)') 'closehisto plot 0'
6121 c-----------------------------------------------------------------------
6122 subroutine xEmsI2(iii,kc)
6123 c-----------------------------------------------------------------------
6124 c plot quanities vs iter
6125 c plot 1: <x> for Pomeron vs iter
6126 c plot 2: <x> for projectile vs iter
6127 c plot 3: <x> for target vs iter
6130 c kc: iteration step
6131 c omega: config probability
6132 c-----------------------------------------------------------------------
6135 include 'epos.incems'
6138 common/cmc1/xp(0:nbim),xt(0:nbim),x(0:nbim),o(0:nbim)
6141 double precision xp,xt,x,xpo,xpj,xtg
6142 common/cemsi2/xpo,xpj,xtg
6151 if(nprmx(k).gt.0)then
6153 if(idpr(n,k).gt.0.and.ivpr(n,k).gt.0)then
6160 if(npom.gt.0)xpo=xpo/npom
6165 if(xpp(i).lt.0.999)then
6166 xpj=xpj+xpp(i)!*xmp(i)
6170 if(npk.gt.0)xpj=xpj/dble(npk)
6175 if(xmt(j).lt.0.999)then
6176 xtg=xtg+xmt(j)!*xpt(j)
6180 if(ntk.gt.0)xtg=xtg/dble(ntk)
6186 elseif(iii.eq.2)then
6191 write(ifhi,'(a)') '!##################################'
6192 write(ifhi,'(a,i3)') '! average x Pom for event ',nrevt+1
6193 write(ifhi,'(a)') '!##################################'
6194 write(ifhi,'(a,i1)') 'openhisto name avxPom-',nrevt+1
6195 write(ifhi,'(a)') 'htyp lin'
6196 write(ifhi,'(a)') 'xmod lin ymod lin'
6197 write(ifhi,'(a,2e11.3)')'xrange',x1,x2
6198 write(ifhi,'(a)') 'text 0 0 "xaxis iteration"'
6199 write(ifhi,'(a)') 'text 0 0 "yaxis average x Pomeron"'
6200 write(ifhi,'(a)') 'array 2'
6202 write(ifhi,'(2e11.3)')float(k),x(k)
6204 write(ifhi,'(a)') ' endarray'
6205 write(ifhi,'(a)') 'closehisto plot 0'
6207 write(ifhi,'(a)') '!##################################'
6208 write(ifhi,'(a,i3)') '! average x proj for event ',nrevt+1
6209 write(ifhi,'(a)') '!##################################'
6210 write(ifhi,'(a,i1)') 'openhisto name avxProj-',nrevt+1
6211 write(ifhi,'(a)') 'htyp lin'
6212 write(ifhi,'(a)') 'xmod lin ymod lin'
6213 write(ifhi,'(a,2e11.3)')'xrange',x1,x2
6214 write(ifhi,'(a)') 'text 0 0 "xaxis iteration"'
6215 write(ifhi,'(a)') 'text 0 0 "yaxis average x proj"'
6216 write(ifhi,'(a)') 'array 2'
6218 write(ifhi,'(2e11.3)')float(k),xp(k)
6220 write(ifhi,'(a)') ' endarray'
6221 write(ifhi,'(a)') 'closehisto plot 0'
6223 write(ifhi,'(a)') '!##################################'
6224 write(ifhi,'(a,i3)') '! average x targ for event ',nrevt+1
6225 write(ifhi,'(a)') '!##################################'
6226 write(ifhi,'(a,i1)') 'openhisto name avxTarg-',nrevt+1
6227 write(ifhi,'(a)') 'htyp lin'
6228 write(ifhi,'(a)') 'xmod lin ymod lin'
6229 write(ifhi,'(a,2e11.3)')'xrange',x1,x2
6230 write(ifhi,'(a)') 'text 0 0 "xaxis iteration"'
6231 write(ifhi,'(a)') 'text 0 0 "yaxis average x targ"'
6232 write(ifhi,'(a)') 'array 2'
6234 write(ifhi,'(2e11.3)')float(k),xt(k)
6236 write(ifhi,'(a)') ' endarray'
6237 write(ifhi,'(a)') 'closehisto plot 0'
6243 c-----------------------------------------------------------------------
6244 subroutine xEmsRx(iii,id,xp,xm)
6245 c-----------------------------------------------------------------------
6246 c plot x+, x-, x, y distribution of remnants
6247 c-----------------------------------------------------------------------
6251 parameter(nbix=50,nbiy=50,nid=2)
6252 common/cxp/nxp(nid),nxm(nid),nx(nid),ny(nid)
6253 *,wxp(nbix,nid),wxm(nbix,nid),wx(nbix,nid),wy(nbiy,nid)
6254 *,xpu,xpo,xmu,xmo,xu,xo,yu,yo,dy
6256 if(iemsrx.eq.0)call utstop('ERROR in XemsRx: iemsrx = 0&')
6284 elseif(iii.eq.1)then
6288 i=1+int(alog(xp/xpu)/alog(xpo/xpu)*nbix)
6291 wxp(i,id)=wxp(i,id)+1
6296 i=1+int(alog(xm/xmu)/alog(xmo/xmu)*nbix)
6299 wxm(i,id)=wxm(i,id)+1
6305 i=1+int(alog(x/xu)/alog(xo/xu)*nbix)
6323 elseif(iii.eq.2)then
6326 if(j.eq.1)iclrem=iclpro
6327 if(j.eq.2)iclrem=icltar
6328 write(ifhi,'(a)') '!----------------------------------'
6329 write(ifhi,'(a)') '! remnant xp distribution '
6330 write(ifhi,'(a)') '!----------------------------------'
6331 write(ifhi,'(a,i1)') 'openhisto name xpRemnant-',j
6332 write(ifhi,'(a)') 'htyp lin'
6333 write(ifhi,'(a)') 'xmod log ymod log'
6334 write(ifhi,'(a,2e11.3)')'xrange',xpu,xpo
6335 write(ifhi,'(a)') 'text 0 0 "xaxis remnant x+"'
6336 write(ifhi,'(a)') 'text 0 0 "yaxis P(x+)"'
6337 write(ifhi,'(a)') 'array 2'
6339 x=xpu*(xpo/xpu)**((i-0.5)/nbix)
6340 dx=xpu*(xpo/xpu)**(1.*i/nbix)*(1.-(xpo/xpu)**(-1./nbix))
6341 if(nxp(j).ne.0)write(ifhi,'(2e11.3)')x,wxp(i,j)/dx/nxp(j)
6342 if(nxp(j).eq.0)write(ifhi,'(2e11.3)')x,0.
6344 write(ifhi,'(a)') ' endarray'
6345 write(ifhi,'(a)') 'closehisto plot 0-'
6346 write(ifhi,'(a)') 'openhisto'
6347 write(ifhi,'(a)') 'htyp lin'
6348 write(ifhi,'(a)') 'array 2'
6350 x=xu*(xo/xu)**((i-0.5)/nbix)
6351 write(ifhi,'(2e11.3)')x,x**alplea(iclrem)*(1+alplea(iclrem))
6353 write(ifhi,'(a)') ' endarray'
6354 write(ifhi,'(a)') 'closehisto plot 0'
6356 write(ifhi,'(a)') '!----------------------------------'
6357 write(ifhi,'(a)') '! remnant xm distribution '
6358 write(ifhi,'(a)') '!----------------------------------'
6359 write(ifhi,'(a,i1)') 'openhisto name xmRemnant-',j
6360 write(ifhi,'(a)') 'htyp lin'
6361 write(ifhi,'(a)') 'xmod log ymod log'
6362 write(ifhi,'(a,2e11.3)')'xrange',xmu,xmo
6363 write(ifhi,'(a)') 'text 0 0 "xaxis remnant x-"'
6364 write(ifhi,'(a)') 'text 0 0 "yaxis P(x-)"'
6365 write(ifhi,'(a)') 'array 2'
6367 x=xmu*(xmo/xmu)**((i-0.5)/nbix)
6368 dx=xmu*(xmo/xmu)**(1.*i/nbix)*(1.-(xmo/xmu)**(-1./nbix))
6369 if(nxm(j).ne.0)write(ifhi,'(2e11.3)')x,wxm(i,j)/dx/nxm(j)
6370 if(nxm(j).eq.0)write(ifhi,'(2e11.3)')x,0.
6372 write(ifhi,'(a)') ' endarray'
6373 write(ifhi,'(a)') 'closehisto plot 0'
6375 write(ifhi,'(a)') '!----------------------------------'
6376 write(ifhi,'(a)') '! remnant x distribution '
6377 write(ifhi,'(a)') '!----------------------------------'
6378 write(ifhi,'(a,i1)') 'openhisto name xRemnant-',j
6379 write(ifhi,'(a)') 'htyp lin'
6380 write(ifhi,'(a)') 'xmod log ymod log'
6381 write(ifhi,'(a,2e11.3)')'xrange',xu,xo
6382 write(ifhi,'(a)') 'text 0 0 "xaxis remnant x"'
6383 write(ifhi,'(a)') 'text 0 0 "yaxis P(x)"'
6384 write(ifhi,'(a)') 'array 2'
6386 x=xu*(xo/xu)**((i-0.5)/nbix)
6387 dx=xu*(xo/xu)**(1.*i/nbix)*(1.-(xo/xu)**(-1./nbix))
6388 if(nx(j).ne.0)write(ifhi,'(2e11.3)')x,wx(i,j)/dx/nx(j)
6389 if(nx(j).eq.0)write(ifhi,'(2e11.3)')x,0.
6391 write(ifhi,'(a)') ' endarray'
6392 write(ifhi,'(a)') 'closehisto plot 0'
6394 write(ifhi,'(a)') '!----------------------------------'
6395 write(ifhi,'(a)') '! remnant y distribution '
6396 write(ifhi,'(a)') '!----------------------------------'
6397 write(ifhi,'(a,i1)') 'openhisto name yRemnant-',j
6398 write(ifhi,'(a)') 'htyp lin'
6399 write(ifhi,'(a)') 'xmod lin ymod log'
6400 write(ifhi,'(a,2e11.3)')'xrange',yu,yo
6401 write(ifhi,'(a)') 'text 0 0 "xaxis remnant y"'
6402 write(ifhi,'(a)') 'text 0 0 "yaxis P(y)"'
6403 write(ifhi,'(a)') 'array 2'
6406 if(ny(j).ne.0)write(ifhi,'(2e11.3)')y,wy(i,j)/dy/ny(j)
6407 if(ny(j).eq.0)write(ifhi,'(2e11.3)')y,0.
6409 write(ifhi,'(a)') ' endarray'
6410 write(ifhi,'(a)') 'closehisto plot 0'
6419 c-----------------------------------------------------------------------
6420 subroutine xEmsPm(iii,ko,nmc)
6421 c-----------------------------------------------------------------------
6422 c m (pomeron number) distribution for different b-bins.
6424 c iii: modus (0,1,2)
6425 c ko: pair number (1 - AB)
6426 c nmc: number of pomerons
6427 c-----------------------------------------------------------------------
6429 include 'epos.incems'
6430 common/geom/rmproj,rmtarg,bmax,bkmx
6431 parameter(nbin=npommx)
6433 common/cn/wn(0:nbin,nbib),wnmc(0:nbin,nbib),npmx(nbib),nn(nbib)
6435 common/cb1/db,b1,b2,bb(nbib),nbibx
6436 double precision plc,s,om1intbc
6439 common/cemspm/sumb(nbib)
6441 if(iemspm.eq.0)call utstop('ERROR in XemsPm: iemspm = 0&')
6459 elseif(iii.eq.1)then
6461 k=int((bk(ko)-b1)/db)+1
6462 if(k.gt.nbibx)k=nbibx
6464 if(nmc.gt.nbin)return
6467 wnmc(nmc,k)=wnmc(nmc,k)+1
6468 sumb(k)=sumb(k)+bk(ko)
6471 elseif(iii.eq.2)then
6476 if(maproj.eq.1.and.matarg.eq.1.and.bmaxim.eq.0.)bb(k)=b1
6477 om1i=sngl(om1intbc(bb(k)))
6482 wn(i,k)=wn(i-1,k)*om1i/i
6484 if(wn(i,k).gt.0.000001*(1.-exp(-om1i)))npmx(k)=i
6487 write(ifhi,'(a)') '!##################################'
6488 write(ifhi,'(a)') '! distr of Pomeron number vs b'
6489 write(ifhi,'(a)') '!##################################'
6490 write(ce,'(f8.2)')sngl(plc)
6491 write(cb,'(f4.2)')bb(k)
6493 write(ifhi,'(a,i1)') 'openhisto name mPom-',k
6494 write(ifhi,'(a)') 'htyp lru'
6495 write(ifhi,'(a)') 'xmod lin ymod log'
6496 write(ifhi,'(a,2e11.3)')'xrange',0.,float(npmx(k))
6497 write(ifhi,'(a)') 'text 0 0 "xaxis number m of Pomerons"'
6498 write(ifhi,'(a)') 'text 0 0 "yaxis prob(m)"'
6500 *write(ifhi,'(a,a)') 'text 0.5 0.90 "E ='//ce//'"'
6501 write(ifhi,'(a,a)') 'text 0.5 0.80 "b ='//cb//'"'
6502 write(ifhi,'(a)') 'array 2'
6504 write(ifhi,'(2e11.3)')float(i),wnmc(i,k)/max(1,nn(k))
6506 write(ifhi,'(a)') ' endarray'
6507 write(ifhi,'(a)') 'closehisto plot 0-'
6510 write(ifhi,'(a)') '!##################################'
6511 write(ifhi,'(a)') '! distr of Pomeron number vs b'
6512 write(ifhi,'(a)') '! traditional approach'
6513 write(ifhi,'(a)') '!##################################'
6514 write(ifhi,'(a,i1)') 'openhisto name mPomTradi-',k
6515 write(ifhi,'(a)') 'htyp lba'
6516 write(ifhi,'(a)') 'xmod lin ymod log'
6517 write(ifhi,'(a,2e11.3)')'xrange',0.,float(npmx(k))
6518 write(ifhi,'(a)') 'array 2'
6520 write(ifhi,'(2e11.3)')float(i),wn(i,k)
6522 write(ifhi,'(a)') ' endarray'
6523 write(ifhi,'(a)') 'closehisto plot 0'
6533 c-----------------------------------------------------------------------
6534 subroutine xEmsB(iii,jjj,ko)
6535 c-----------------------------------------------------------------------
6536 c b distribution at different stages
6538 c iii: modus (0,1,2)
6539 c jjj: stage or type of interaction
6540 c just after Metropolis:
6543 c after defining diffraction:
6548 c ko: pair number (1 - AB)
6549 c-----------------------------------------------------------------------
6551 include 'epos.incems'
6552 include 'epos.incsem'
6555 common/cxemsb1/w(0:njjj,nbib),nn(njjj)
6556 common/cxemsb2/db,b1,b2
6557 common/cxemsb3/njjj1
6558 double precision PhiExact,om1intbi,PhiExpo!,PhiUnit
6559 common/geom/rmproj,rmtarg,bmax,bkmx
6560 dimension uua2(nbib),uuo2(nbib),uu3(nbib)
6562 if(iemsb.eq.0)call utstop('ERROR in XemsB: iemsB = 0&')
6576 elseif(iii.eq.1)then
6581 k=int((bk(ko)-b1)/db)+1
6588 elseif(iii.eq.2)then
6590 if(njjj1.ne.1)call utstop
6591 &('xEmsB must be called also with jjj=1&')
6597 y=w(1,k)/nn(1)/(pi*((x+0.5*db)**2-(x-0.5*db)**2))
6606 write(ifhi,'(a)') '!##################################'
6607 write(ifhi,'(a)') '! b distr exact theory '
6608 write(ifhi,'(a)') '!##################################'
6609 if(j.ge.2.and.j.le.6)then
6610 write(ifhi,'(a,i1,a)') 'openhisto name b',j,'Exact'
6611 write(ifhi,'(a)') 'htyp lba xmod lin ymod lin'
6612 write(ifhi,'(a)') 'text 0 0 "xaxis impact parameter b"'
6613 write(ifhi,'(a)') 'text 0 0 "yaxis P(b)"'
6614 write(ifhi,'(a)') 'array 2'
6618 uuo2(k)=sngl(PhiExpo(1.,1.d0,1.d0,engy**2,b))
6619 uua2(k)=min(uuo2(k),max(0.,
6620 & sngl(PhiExact(1.,1.d0,1.d0,engy**2,b))))
6621 uu3(k)=sngl(min(50d0,exp(om1intbi(b,2)/dble(r2hads(iclpro)
6622 & +r2hads(icltar)))))
6624 if(j.eq.2)y=(1.-uua2(k))
6626 if(j.eq.4)y=(1.-uua2(k)*uu3(k))
6627 if(j.eq.5)y=uua2(k)*(uu3(k)-1.)
6628 if(j.eq.6)y=(1.-uua2(k))
6629 write(ifhi,'(2e11.3)')b,y
6631 write(ifhi,'(a)') ' endarray'
6632 write(ifhi,'(a)') 'closehisto plot 0-'
6634 write(ifhi,'(a)') '!##################################'
6635 write(ifhi,'(a)') '! b distr unitarized theory '
6636 write(ifhi,'(a)') '!##################################'
6637 write(ifhi,'(a,i1,a)') 'openhisto name b',j,'Unit'
6638 write(ifhi,'(a)') 'htyp lbf xmod lin ymod lin'
6639 write(ifhi,'(a)') 'text 0 0 "xaxis impact parameter b"'
6640 write(ifhi,'(a)') 'text 0 0 "yaxis P(b)"'
6641 write(ifhi,'(a)') 'array 2'
6645 if(j.eq.2)y=(1.-uuo2(k))
6647 if(j.eq.4)y=(1.-uuo2(k)*uu3(k))
6648 if(j.eq.5)y=uuo2(k)*(uu3(k)-1.)
6649 if(j.eq.6)y=(1.-uuo2(k))
6650 write(ifhi,'(2e11.3)')b,y
6652 write(ifhi,'(a)') ' endarray'
6653 write(ifhi,'(a)') 'closehisto plot 0-'
6654 write(ifhi,'(a)') '!##################################'
6655 write(ifhi,'(a)') '! b distr for cross section '
6656 write(ifhi,'(a)') '!##################################'
6657 write(ifhi,'(a,i1,a)') 'openhisto name b',j,'Unit'
6658 write(ifhi,'(a)') 'htyp lge xmod lin ymod lin'
6659 write(ifhi,'(a)') 'text 0 0 "xaxis impact parameter b"'
6660 write(ifhi,'(a)') 'text 0 0 "yaxis P(b)"'
6661 write(ifhi,'(a)') 'array 2'
6665 if(j.eq.2)y=(1.-(uuo2(k)+uua2(k))*0.5)
6666 if(j.eq.3)y=(uuo2(k)+uua2(k))*0.5
6667 if(j.eq.4)y=(1.-(uuo2(k)+uua2(k))*0.5*uu3(k))
6668 if(j.eq.5)y=(uuo2(k)+uua2(k))*0.5*(uu3(k)-1.)
6669 if(j.eq.6)y=(1.-(uuo2(k)+uua2(k))*0.5)
6670 write(ifhi,'(2e11.3)')b,y
6672 write(ifhi,'(a)') ' endarray'
6673 write(ifhi,'(a)') 'closehisto plot 0-'
6674 write(ifhi,'(a)') '!##################################'
6675 write(ifhi,'(a)') '! b distribution simulation'
6676 write(ifhi,'(a)') '!##################################'
6677 write(ifhi,'(a,i1,a)') 'openhisto name b',j,'Simu'
6678 write(ifhi,'(a)') 'htyp lrf xmod lin ymod lin'
6679 write(ifhi,'(a,2e11.3)')'xrange',0.0,b2
6680 write(ifhi,'(a,2e11.3)')'yrange',0.,ymax
6681 write(ifhi,'(a)') 'text 0 0 "xaxis impact parameter b"'
6682 write(ifhi,'(a)') 'text 0 0 "yaxis P(b)"'
6683 if(j.eq.1)write(ifhi,'(a)')'text 0.1 0.35 "after Metropolis"'
6684 if(j.eq.1)write(ifhi,'(a)')'text 0.2 0.20 "all "'
6685 if(j.eq.2)write(ifhi,'(a)')'text 0.3 0.85 "after Metropolis"'
6686 if(j.eq.2)write(ifhi,'(a)')'text 0.5 0.70 "interaction "'
6687 if(j.eq.3)write(ifhi,'(a)')'text 0.3 0.85 "nothing"'
6688 if(j.eq.4)write(ifhi,'(a)')'text 0.3 0.85 "cut"'
6689 if(j.eq.5)write(ifhi,'(a)')'text 0.3 0.85 "diffr"'
6690 if(j.eq.6)write(ifhi,'(a)')'text 0.3 0.85 "cut + diffr"'
6691 write(ifhi,'(a)') 'array 2'
6694 if(j.eq.1)y=fk*w(j,k)/nn(1)/(pi*((x+0.5*db)**2-(x-0.5*db)**2))
6696 if(j.ne.1.and.w(1,k).ne.0.)y=w(j,k)/w(1,k)
6697 if(nn(j).gt.0)write(ifhi,'(2e11.3)')x,y
6699 write(ifhi,'(a)') ' endarray'
6700 write(ifhi,'(a)') 'closehisto plot 0'
6711 c-----------------------------------------------------------------------
6712 subroutine xEmsBg(iii,jjj,ko)
6713 c-----------------------------------------------------------------------
6714 c b distribution at different stages for different group
6716 c iii: modus (0,1,2,3)
6717 c jjj: group of interaction (1,2 ... ,7)
6718 c ko: pair number (1 - AB)
6719 c-----------------------------------------------------------------------
6721 include 'epos.incems'
6724 common/cxemsb4/wg(-1:njjj,nbib),nng(nbib),uug(nbib),kollx
6725 common/cxemsb5/dbg,b1g,b2g
6726 common/cxemsb6/njjj0
6727 double precision seedp,PhiExpo!,PhiExact
6728 common/geom/rmproj,rmtarg,bmax,bkmx
6730 if(iemsbg.eq.0)call utstop('ERROR in XemsBg: iemsbg = 0&')
6743 elseif(iii.eq.1)then
6748 k=int((bk(ko)-b1g)/dbg)+1
6751 if(jjj.eq.-1.or.jjj.eq.0)then
6752 wg(jjj,k)=wg(jjj,k)+1
6754 wg(jjj,k)=wg(jjj,k)+1
6759 elseif(iii.eq.3)then
6763 om1i=sngl(om1intc(k))
6764 if(rangen().lt.1.-exp(-om1i))then
6765 c om1i=sngl(PhiExpo(1.,1.d0,1.d0,engy*engy,bk(k)))
6766 c if(rangen().lt.1.-om1i)then
6772 elseif(iii.eq.2)then
6774 if(njjj0.ne.1)call utstop
6775 &('xEmsBg must be called also with jjj=0&')
6781 if(matarg+maproj.gt.2)then
6786 wtot=wtot/float(kollx)
6791 write(ifhi,'(a)') '!##################################'
6792 write(ifhi,'(a)') '! b distribution simulation'
6793 write(ifhi,'(a)') '!##################################'
6794 write(ifhi,'(a,i1,a)') 'openhisto name bg',j,'Simu'
6795 write(ifhi,'(a)') 'htyp lin xmod lin ymod lin'
6796 write(ifhi,'(a,2e11.3)')'xrange',0.,b2g
6797 write(ifhi,'(a,2e11.3)')'yrange',0.,ymax
6798 write(ifhi,'(a)') 'text 0 0 "xaxis impact parameter b"'
6799 write(ifhi,'(a)') 'text 0 0 "yaxis P(b)"'
6801 &write(ifhi,'(a,f7.4,a)') 'text 0.5 0.8 "alpha=',1./wtot,'"'
6802 write(ifhi,'(a)') 'array 2'
6806 if(nng(k).ne.0.and.wg(0,k).ne.0)
6807 & y=wg(j,k)/float(nng(k))*wg(-1,k)/wg(0,k)!/wtot
6808 c if(wg(0,k).ne.0..and.nng(k).ne.0)y=wg(j,k)/nng(k)*wg(-1,k)/wg(0,k)
6809 c!???????????? better normalization ? probability to have an interaction
6810 c in epos compared to eikonal probability, instead of normalized by the
6811 c probability of a collision for a pair (the number collision/number
6814 write(ifhi,'(2e11.3)')b,y
6816 write(ifhi,'(a)') ' endarray'
6817 write(ifhi,'(a)') 'closehisto plot 0-'
6819 write(ifhi,'(a)') '!##################################'
6820 write(ifhi,'(a)') '! b distr tot simul theory '
6821 write(ifhi,'(a)') '!##################################'
6822 write(ifhi,'(a)') 'openhisto name btotSimu'
6823 write(ifhi,'(a)') 'htyp pfc xmod lin ymod lin'
6824 write(ifhi,'(a)') 'text 0 0 "xaxis impact parameter b"'
6825 write(ifhi,'(a)') 'text 0 0 "yaxis P(b)"'
6826 write(ifhi,'(a)') 'array 2'
6829 write(ifhi,'(2e11.3)')b,uug(k)
6831 write(ifhi,'(a)') ' endarray'
6832 write(ifhi,'(a)') 'closehisto plot 0-'
6833 write(ifhi,'(a)') '!##################################'
6834 write(ifhi,'(a)') '! b distr unitarized theory '
6835 write(ifhi,'(a)') '!##################################'
6836 write(ifhi,'(a,i1,a)') 'openhisto name bg',j,'Unit'
6837 write(ifhi,'(a)') 'htyp lba xmod lin ymod lin'
6838 write(ifhi,'(a)') 'text 0 0 "xaxis impact parameter b"'
6839 write(ifhi,'(a)') 'text 0 0 "yaxis P(b)"'
6840 write(ifhi,'(a)') 'array 2'
6843 c a1=PhiExact(1.,1.d0,1.d0,engy**2,b)
6844 a1=sngl(PhiExpo(1.,1.d0,1.d0,engy**2,b))
6846 write(ifhi,'(2e11.3)')b,y
6848 write(ifhi,'(a)') ' endarray'
6849 write(ifhi,'(a)') 'closehisto plot 0'
6858 c-----------------------------------------------------------------------
6859 subroutine xEmsPx(iii,xmc,ymc,npos)
6860 c-----------------------------------------------------------------------
6861 c plot x-distribution and y-distribution of Pomerons
6862 c-----------------------------------------------------------------------
6865 include 'epos.incems'
6866 common/geom/rmproj,rmtarg,bmax,bkmx
6868 parameter(nbix=30,nbib=51)
6869 common/cx/x(2,nbix),dx(2,nbix),wxmc(2,nbix),wxmcI(2,nbix)
6870 * ,xl(2,nbix),dxl(2,nbix),wxp(2,nbix),wxm(2,nbix),wxpI(2,nbix)
6871 *,wxmI(2,nbix),wxpY(2,nbix),wxmY(2,nbix),wxmcY(2,nbix)
6873 common/cy/y(nbiy),wymc(nbiy),wymcY(nbiy),wymcI(nbiy),nyp,nym
6874 double precision PomIncXExact,PomIncPExact,PomIncMExact,dcel
6875 double precision PomIncXIExact,PomIncPIExact,PomIncMIExact
6877 common/cemspx/xu,xo,yu,yo,dy,xlu,xlo,bb,nn,db,mm,nm,nt
6878 character mod*5, imod*5, txtxm*6
6882 if(iemspx.eq.0)call utstop('ERROR in XemsPx: iemspx = 0&')
6894 x(1,i)=xu*(xo/xu)**((i-0.5)/nbix)
6895 x(2,i)=xu+(xo-xu)*((i-0.5)/nbix)
6896 dx(1,i)=xu*(xo/xu)**(1.*i/nbix)*(1.-(xo/xu)**(-1./nbix))
6897 dx(2,i)=(xo-xu)/nbix
6906 xl(1,i)=xlu*(xlo/xlu)**((i-0.5)/nbix)
6907 xl(2,i)=xlu+(xlo-xlu)*((i-0.5)/nbix)
6908 dxl(1,i)=xlu*(xlo/xlu)**(1.*i/nbix)*(1.-(xlo/xlu)**(-1./nbix))
6909 dxl(2,i)=(xlo-xlu)/nbix
6924 y(i)=yu+dy/2.+float(i-1)*dy
6933 db=bkmx*2./float(nbib-1)
6935 elseif(iii.eq.1)then
6937 xp=sqrt(xmc)*exp(ymc)
6938 xm=sqrt(xmc)*exp(-ymc)
6942 i=1+int(alog(xmc/xu)/alog(xo/xu)*nbix)
6945 wxmc(1,i)=wxmc(1,i)+1.
6946 if(npos.eq.1) wxmcI(1,i)=wxmcI(1,i)+1.
6947 if(npos.eq.nposi)wxmcY(1,i)=wxmcY(1,i)+1.
6949 i=1+int((xmc-xu)/(xo-xu)*nbix)
6952 wxmc(2,i)=wxmc(2,i)+1.
6953 if(npos.eq.1) wxmcI(2,i)=wxmcI(2,i)+1.
6954 if(npos.eq.nposi)wxmcY(2,i)=wxmcY(2,i)+1.
6958 i=1+int(alog(xp/xlu)/alog(xlo/xlu)*nbix)
6961 wxp(1,i)=wxp(1,i)+1.
6962 if(npos.eq.1) wxpI(1,i)=wxpI(1,i)+1.
6963 if(npos.eq.nposi)wxpY(1,i)=wxpY(1,i)+1.
6965 i=1+int((xp-xlu)/(xlo-xlu)*nbix)
6968 wxp(2,i)=wxp(2,i)+1.
6969 if(npos.eq.1) wxpI(2,i)=wxpI(2,i)+1.
6970 if(npos.eq.nposi)wxpY(2,i)=wxpY(2,i)+1.
6974 i=1+int(alog(xm/xlu)/alog(xlo/xlu)*nbix)
6977 wxm(1,i)=wxm(1,i)+1.
6978 if(npos.eq.1) wxmI(1,i)=wxmI(1,i)+1.
6979 if(npos.eq.nposi)wxmY(1,i)=wxmY(1,i)+1.
6981 i=1+int((xm-xlu)/(xlo-xlu)*nbix)
6984 wxm(2,i)=wxm(2,i)+1.
6985 if(npos.eq.1) wxmI(2,i)=wxmI(2,i)+1.
6986 if(npos.eq.nposi)wxmY(2,i)=wxmY(2,i)+1.
6990 i=int((ymc-yu)/dy)+1
6994 if(npos.eq.1) wymcI(i)=wymcI(i)+1
6995 if(npos.eq.nposi)wymcY(i)=wymcY(i)+1
6996 if(ymc.gt.0)nyp=nyp+1
6997 if(ymc.lt.0)nym=nym+1
6999 elseif(iii.eq.2)then
7001 if(maproj.eq.1.and.matarg.eq.1.and.bminim.eq.bmaxim)then
7004 ff=float(nrevt)/float(ntevt)
7006 elseif(maproj.eq.1.and.matarg.eq.1)then
7010 elseif(bminim.lt.0.001.and.bmaxim.gt.20)then
7012 area=pi*(rmproj+rmtarg)**2
7013 ff=area*float(nrevt)/float(ntevt)/(maproj*matarg)/sigine*10
7016 write(ifmt,*)'xEmsPx ignored'
7025 if(kk.eq.1)mod=' log '
7026 if(kk.eq.2)mod=' lin '
7028 write(ifhi,'(a)') '!----------------------------------'
7029 write(ifhi,'(a)') '! Pomeron x distribution '//mod
7030 write(ifhi,'(a)') '!----------------------------------'
7032 write(ifhi,'(a)') 'openhisto name xPomSimuL'//mod(3:4)
7033 write(ifhi,'(a)') 'htyp lru xmod'//mod//'ymod log'
7034 write(ifhi,'(a,2e11.3)')'xrange',xu,xo
7035 write(ifhi,'(a)') 'text 0 0 "xaxis x?PE!"'
7036 write(ifhi,'(a)') 'text 0 0 "yaxis'//imod//'?Pom! / dx?PE!"'
7037 if(kk.eq.1)write(ifhi,'(a,f5.2,a)')'text 0.1 0.3 "f=',ff,'"'
7038 if(kk.eq.2)write(ifhi,'(a,f5.2,a)')'text 0.1 0.1 "f=',ff,'"'
7039 write(ifhi,'(a)') 'array 2'
7043 z=ff*wxmc(kk,i)/dx(kk,i)/nrevt
7045 write(ifhi,'(2e11.3)')u,z
7047 write(ifhi,'(a)') ' endarray'
7048 write(ifhi,'(a)') 'closehisto plot 0-'
7050 write(ifhi,'(a)') 'openhisto name xPomUnitL'//mod(3:4)
7051 write(ifhi,'(a)') 'htyp lba xmod'//mod//'ymod log'
7052 write(ifhi,'(a,2e11.3)')'xrange',xu,xo
7053 write(ifhi,'(a)') 'text 0 0 "xaxis x?PE!"'
7054 write(ifhi,'(a)') 'text 0 0 "yaxis'//imod//'?Pom! / dx?PE!"'
7055 write(ifhi,'(a)') 'array 2'
7059 if(mmmm.eq.1)z=PomIncXExact(dble(u),bb)
7060 if(mmmm.eq.2)z=PomIncXIExact(dble(u))/sigine*10
7061 if(mmmm.eq.3)z=PomIncXIExact(dble(u))/sigine*10
7063 write(ifhi,'(2e11.3)')u,z
7065 write(ifhi,'(a)') ' endarray'
7066 write(ifhi,'(a,f5.3,a,f5.3,a)')
7067 * 'text .1 .85 "I= ',s1,' (',s2,')"'
7068 write(ifhi,'(a)') 'closehisto plot 0'
7070 write(ifhi,'(a)') '!--------------------------------'
7071 write(ifhi,'(a)') '! Pomeron y distribution '//mod
7072 write(ifhi,'(a)') '!--------------------------------'
7074 write(ifhi,'(a)') 'openhisto name yPomSimuL'//mod(3:4)
7075 write(ifhi,'(a)') 'htyp lru xmod lin ymod'//mod
7076 write(ifhi,'(a,2e11.3)')'xrange',yu,yo
7077 write(ifhi,'(a)') 'text 0 0 "xaxis y?PE!"'
7078 write(ifhi,'(a)') 'text 0 0 "yaxis'//imod//'?Pom!/dy?PE!"'
7079 write(ifhi,'(a,f5.2,a)')'text 0.1 0.7 "f=',ff,'"'
7080 write(ifhi,'(a)') 'array 2'
7084 z=ff*wymc(i)/dy/nrevt
7086 write(ifhi,'(2e11.3)')u,z
7088 write(ifhi,'(a)') ' endarray'
7089 write(ifhi,'(a)') 'closehisto plot 0'
7091 write(ifhi,'(a)') '!----------------------------------'
7092 write(ifhi,'(a)') '! Pomeron x+ distribution '//mod
7093 write(ifhi,'(a)') '!----------------------------------'
7095 write(ifhi,'(a)') 'openhisto name xpPomSimuL'//mod(3:4)
7096 write(ifhi,'(a)') 'htyp lru xmod'//mod//'ymod log'
7097 write(ifhi,'(a,2e11.3)')'xrange',xlu,xlo
7098 write(ifhi,'(a)') 'text 0 0 "xaxis x+?PE!"'
7099 write(ifhi,'(a)') 'text 0 0 "yaxis'//imod//'?Pom! / dx+?PE!"'
7100 if(kk.eq.1)write(ifhi,'(a,f5.2,a)')'text 0.1 0.3 "f=',ff,'"'
7101 if(kk.eq.2)write(ifhi,'(a,f5.2,a)')'text 0.1 0.1 "f=',ff,'"'
7102 write(ifhi,'(a)') 'array 2'
7106 z=ff*wxp(kk,i)/dxl(kk,i)/nrevt
7108 write(ifhi,'(2e11.3)')u,z
7110 write(ifhi,'(a)') ' endarray'
7111 write(ifhi,'(a)') 'closehisto plot 0-'
7113 write(ifhi,'(a)') 'openhisto name xpPomUnitL'//mod(3:4)
7114 write(ifhi,'(a)') 'htyp lba xmod'//mod//'ymod log'
7115 write(ifhi,'(a,2e11.3)')'xrange',xlu,xlo
7116 write(ifhi,'(a)') 'text 0 0 "xaxis x+?PE!"'
7117 write(ifhi,'(a)') 'text 0 0 "yaxis'//imod//'?Pom! / dx+?PE!"'
7118 write(ifhi,'(a)') 'array 2'
7122 if(mmmm.eq.1)z=PomIncPExact(dble(u),bb)
7123 if(mmmm.eq.2)z=PomIncPIExact(dble(u))/sigine*10
7124 if(mmmm.eq.3)z=PomIncPIExact(dble(u))/sigine*10
7126 write(ifhi,'(2e11.3)')u,z
7128 write(ifhi,'(a)') ' endarray'
7129 write(ifhi,'(a,f5.3,a,f5.3,a)')
7130 * 'text .1 .85 "I= ',s1,' (',s2,')"'
7131 write(ifhi,'(a)') 'closehisto plot 0'
7133 write(ifhi,'(a)') '!----------------------------------'
7134 write(ifhi,'(a)') '! x-?PE! distribution '//mod
7135 write(ifhi,'(a)') '!----------------------------------'
7137 write(ifhi,'(a)') 'openhisto name xmPomSimuL'//mod(3:4)
7138 write(ifhi,'(a)') 'htyp lru xmod'//mod//'ymod log'
7139 write(ifhi,'(a,2e11.3)')'xrange',xlu,xlo
7140 write(ifhi,'(a)') 'text 0 0 "xaxis x-?PE!"'
7141 write(ifhi,'(a)') 'text 0 0 "yaxis'//imod//'?Pom! / dx-?PE!"'
7142 if(kk.eq.1)write(ifhi,'(a,f5.2,a)')'text 0.1 0.3 "f=',ff,'"'
7143 if(kk.eq.2)write(ifhi,'(a,f5.2,a)')'text 0.1 0.1 "f=',ff,'"'
7144 write(ifhi,'(a)') 'array 2'
7148 z=ff*wxm(kk,i)/dxl(kk,i)/nrevt
7150 write(ifhi,'(2e11.3)')u,z
7152 write(ifhi,'(a)') ' endarray'
7153 write(ifhi,'(a)') 'closehisto plot 0-'
7155 write(ifhi,'(a)') 'openhisto name xmPomUnitL'//mod(3:4)
7156 write(ifhi,'(a)') 'htyp lba xmod'//mod//'ymod log'
7157 write(ifhi,'(a,2e11.3)')'xrange',xlu,xlo
7158 write(ifhi,'(a)') 'text 0 0 "xaxis x-?PE!"'
7159 write(ifhi,'(a)') 'text 0 0 "yaxis'//imod//'?Pom! / dx-"'
7160 write(ifhi,'(a)') 'array 2'
7164 if(mmmm.eq.1)z=PomIncMExact(dble(u),bb)
7165 if(mmmm.eq.2)z=PomIncMIExact(dble(u))/sigine*10
7166 if(mmmm.eq.3)z=PomIncMIExact(dble(u))/sigine*10
7168 write(ifhi,'(2e11.3)')u,z
7170 write(ifhi,'(a)') ' endarray'
7171 write(ifhi,'(a,f5.3,a,f5.3,a)')
7172 * 'text .1 .85 "I= ',s1,' (',s2,')"'
7173 write(ifhi,'(a)') 'closehisto plot 0'
7175 !................................................................
7177 xm=-1. !xm integration
7183 write(ifhi,'(a)') '!----------------------------------'
7184 write(ifhi,'(a,3i1)') '! ffom11 '//mod,jjb,jj
7185 write(ifhi,'(a)') '!----------------------------------'
7187 write(ifhi,'(a,2i1)')'openhisto name ffom11L'//mod(3:4),jjb,jj+8
7188 write(ifhi,'(a)') 'htyp lin xmod'//mod//'ymod log'
7189 write(ifhi,'(a,2e11.3)')'xrange ',xlu,xlo
7190 write(ifhi,'(a)')'txt "xaxis x+?PE!"'
7191 write(ifhi,'(a)')'txt "yaxis dn?Pom! / dx+?PE! "'
7192 write(ifhi,'(a)')'text 0.05 0.1 "fit and exact, all contrib."'
7193 if(jjb.lt.3)write(ifhi,'(a,f4.1,3a)')
7194 * 'txt "title ffom11 b =',b,' ',txtxm,'"'
7195 if(jjb.ge.3)write(ifhi,'(3a)')
7196 * 'txt "title ffom11 b aver ',txtxm,'"'
7197 write(ifhi,'(a)') 'array 2'
7200 if(jjb.lt.3.and.jj.eq.0)z= ffom11(u,xm,b,-1,-1)
7201 if(jjb.lt.3.and.jj.eq.1)z= ffom11(u,xm,b,0,5)
7202 if(jjb.lt.3.and.jj.eq.2)z= ffom11(u,xm,b,0,4)
7203 if(jjb.eq.3.and.jj.eq.0)z=ffom11a(u,xm,-1,-1)
7204 if(jjb.eq.3.and.jj.eq.1)z=ffom11a(u,xm,0,5)
7205 if(jjb.eq.3.and.jj.eq.2)z=ffom11a(u,xm,0,4)
7206 write(ifhi,'(2e11.3)')u,z
7208 write(ifhi,'(a)') ' endarray'
7209 if(jj.le.1)write(ifhi,'(a)') 'closehisto plot 0-'
7210 if(jj.eq.2)write(ifhi,'(a)') 'closehisto plot 0'
7221 write(ifhi,'(a)') '!----------------------------------'
7222 write(ifhi,'(a,3i1)') '! ffom11 '//mod,jjb,jj
7223 write(ifhi,'(a)') '!----------------------------------'
7225 write(ifhi,'(a,3i1)')'openhisto name om1ffL'//mod(3:4),jjb,jj
7226 if(jj.ne.0)write(ifhi,'(a)') 'htyp lin xmod'//mod//'ymod log'
7227 if(jj.eq.0)write(ifhi,'(a)') 'htyp lro xmod'//mod//'ymod log'
7228 write(ifhi,'(a,2e11.3)')'xrange ',xlu,xlo
7230 write(ifhi,'(a)') 'txt "xaxis x+?PE!"'
7231 write(ifhi,'(a)') 'txt "yaxis dn?Pom! / dx+?PE! "'
7233 write(ifhi,'(a)') 'text 0.1 0.2 "soft sea-sea"'
7234 write(ifhi,'(a)') 'text 0.1 0.1 "val-sea sea-val val-val"'
7236 write(ifhi,'(a)') 'text 0.05 0.8 "soft"'
7237 write(ifhi,'(a)') 'text 0.05 0.7 "diff"'
7238 write(ifhi,'(a)') 'text 0.05 0.6 "sea-sea"'
7239 write(ifhi,'(a)') 'text 0.05 0.5 "val-sea"'
7240 write(ifhi,'(a)') 'text 0.05 0.4 "sea-val"'
7241 write(ifhi,'(a)') 'text 0.05 0.3 "val-val"'
7243 if(jjb.lt.3)write(ifhi,'(a,f4.1,3a)')
7244 * 'txt "title ffom11 b =',b,' ',txtxm,'"'
7245 if(jjb.ge.3)write(ifhi,'(3a)')
7246 * 'txt "title ffom11 b aver ',txtxm,'"'
7248 write(ifhi,'(a)') 'array 2'
7251 if(jjb.lt.3)z= ffom11(u,xm,b,jj,jj)
7252 if(jjb.eq.3)z=ffom11a(u,xm,jj,jj)
7253 write(ifhi,'(2e11.3)')u,z
7255 write(ifhi,'(a)') ' endarray'
7256 if(jjj.ne.6)write(ifhi,'(a)') 'closehisto plot 0-'
7257 if(jjj.eq.6)write(ifhi,'(a)') 'closehisto plot 0'
7269 c-----------------------------------------------------------------------
7270 subroutine xEmsP2(iii,jaa,jex,xpd,xmd,xpb,xmb,pt1,pt2)
7271 c-----------------------------------------------------------------------
7272 c plot x+ distributions of Pomeron ends (PE) (xpd)
7273 c and Pomeron's in Born (IB) partons (xpb),
7274 c and pt dist of Pomeron's out Born (OB) partons
7275 c integrated over x- bins (xmd,xmb)
7278 c iii>=2: make histogram
7279 c (2 - Pomeron end PE, 3 - in Born IB, 4 - out Born OB)
7280 c jaa: type of semihard Pomeron
7281 c 1= sea-sea, 2= val=sea, 3= sea-val, 4= val-val
7283 c jex: emission type
7284 c 1= no emission, 2= proj emis, 3= targ emis, 4= both sides
7286 c-----------------------------------------------------------------------
7289 include 'epos.incsem'
7290 include 'epos.incems'
7291 common/geom/rmproj,rmtarg,bmax,bkmx
7292 parameter(nbixp=25,nbixm=5,nbipt=20)
7293 common/cxb/xlp(2,nbixp),dxlp(2,nbixp)
7294 * ,xlm(2,nbixm),dxlm(2,nbixm)
7295 * ,wxb(2,4,4,nbixp,nbixm)
7296 * ,wxe(2,4,4,nbixp,nbixm)
7297 common/cptb/ptu,pto,ptob(nbipt),wptob(4,4,nbipt)
7298 common/cemspbx/xlub1,xlub2,xlob
7299 ctp060829 character imod*5
7301 if(iemspbx.eq.0)call utstop('ERROR in xEmsP2: iemspbx = 0&')
7309 xlp(1,i)=xlub1*(xlob/xlub1)**((i-0.5)/nbixp)
7310 xlp(2,i)=xlub2+(xlob-xlub2)*((i-0.5)/nbixp)
7311 dxlp(1,i)=xlub1*(xlob/xlub1)**(1.*i/nbixp)
7312 * *(1.-(xlob/xlub1)**(-1./nbixp))
7313 dxlp(2,i)=(xlob-xlub2)/nbixp
7316 xlm(1,i)=xlub1*(xlob/xlub1)**((i-0.5)/nbixm)
7317 xlm(2,i)=xlub2+(xlob-xlub2)*((i-0.5)/nbixm)
7318 dxlm(1,i)=xlub1*(xlob/xlub1)**(1.*i/nbixm)
7319 * *(1.-(xlob/xlub1)**(-1./nbixm))
7320 dxlm(2,i)=(xlob-xlub2)/nbixm
7326 wxb(1,jaai,jexi,i,j)=0.
7327 wxb(2,jaai,jexi,i,j)=0.
7328 wxe(1,jaai,jexi,i,j)=0.
7329 wxe(2,jaai,jexi,i,j)=0.
7337 ptob(i)=ptu+(pto-ptu)*(i-0.5)/nbipt
7340 wptob(jaai,jexi,i)=0
7345 elseif(iii.eq.1)then
7349 if(xp.lt.xlub1)goto2
7350 if(xm.lt.xlub1)goto2
7351 i=1+int(alog(xp/xlub1)/alog(xlob/xlub1)*nbixp)
7354 j=1+int(alog(xm/xlub1)/alog(xlob/xlub1)*nbixm)
7357 wxb(1,jaa,jex,i,j)=wxb(1,jaa,jex,i,j)+1.
7360 if(xp.lt.xlub2)goto12
7361 if(xm.lt.xlub2)goto12
7362 i=1+int((xp-xlub2)/(xlob-xlub2)*nbixp)
7363 if(i.gt.nbixp)goto12
7365 j=1+int((xm-xlub2)/(xlob-xlub2)*nbixm)
7366 if(j.gt.nbixm)goto12
7368 wxb(2,jaa,jex,i,j)=wxb(2,jaa,jex,i,j)+1.
7373 if(xp.lt.xlub1)goto22
7374 if(xm.lt.xlub1)goto22
7375 i=1+int(alog(xp/xlub1)/alog(xlob/xlub1)*nbixp)
7376 if(i.gt.nbixp)goto22
7378 j=1+int(alog(xm/xlub1)/alog(xlob/xlub1)*nbixm)
7379 if(j.gt.nbixm)goto22
7381 wxe(1,jaa,jex,i,j)=wxe(1,jaa,jex,i,j)+1.
7384 if(xp.lt.xlub2)goto32
7385 if(xm.lt.xlub2)goto32
7386 i=1+int((xp-xlub2)/(xlob-xlub2)*nbixp)
7387 if(i.gt.nbixp)goto32
7389 j=1+int((xm-xlub2)/(xlob-xlub2)*nbixm)
7390 if(j.gt.nbixm)goto32
7392 wxe(2,jaa,jex,i,j)=wxe(2,jaa,jex,i,j)+1.
7398 i=1+int((pt-ptu)/(pto-ptu)*nbipt)
7400 if(i.gt.nbipt)goto42
7401 wptob(jaa,jex,i)=wptob(jaa,jex,i)+1
7405 elseif(iii.ge.2)then
7407 if(maproj.eq.1.and.matarg.eq.1.and.bminim.eq.bmaxim)then
7410 ff=float(nrevt)/float(ntevt)
7411 ctp060829 imod=' dn'
7412 elseif(maproj.eq.1.and.matarg.eq.1)then
7415 ctp060829 imod=' dn'
7416 elseif(bminim.lt.0.001.and.bmaxim.gt.20)then
7418 area=pi*(rmproj+rmtarg)**2
7419 ff=area*float(nrevt)/float(ntevt)/(maproj*matarg)/sigine*10
7420 ctp060829 imod=' dn'
7422 write(ifmt,*)'xEmsP2 ignored'
7426 j1=1 !nint(xpar1) !first xminus bin
7427 j2=5 !nint(xpar2) !last xminus bin
7429 kkk=2 !nint(xpar3) !1 (log binning) 2 (lin binning)
7431 ctp060829 xmi1=xlub1*(xlob/xlub1)**((j1-1.)/nbixm)
7432 ctp060829 xmi2=xlub1*(xlob/xlub1)**((j2-0.)/nbixm)
7434 elseif(kkk.eq.2)then
7435 ctp060829 xmi1=xlub2+(xlob-xlub2)*((j1-1.)/nbixm)
7436 ctp060829 xmi2=xlub2+(xlob-xlub2)*((j2-0.)/nbixm)
7456 elseif(jex.eq.2)then
7459 elseif(jex.eq.3)then
7462 elseif(jex.eq.4)then
7465 elseif(jex.eq.5)then
7472 write(ifhi,'(a)') '!----------------------------------'
7473 write(ifhi,'(a,3i1)') '! PE ',jaa,jex
7474 write(ifhi,'(a)') '!----------------------------------'
7476 sum=ffom12aii(jaa,je1,je2)
7477 write(ifhi,'(a,2i1)')'openhisto name ffom12a',jaa,jex
7478 write(ifhi,'(a)')'htyp lin xmod lin ymod log'
7479 write(ifhi,'(a,2e11.3)')'xrange ',xlub,xlob
7480 write(ifhi,'(a)') 'txt "xaxis x+?PE!"'
7481 write(ifhi,'(a)') 'txt "yaxis dn?semi! / dx+?PE! "'
7482 write(ifhi,'(a,2i1,a)')'txt "title ffom12a + MC (',jaa,jex,')"'
7483 write(ifhi,'(a)') 'array 2'
7486 z=ffom12ai(u,jaa1,jaa2,je1,je2)
7487 write(ifhi,'(2e11.3)')u,z
7489 write(ifhi,'(a)') ' endarray'
7491 write(ifhi,'(a)') 'closehisto plot 0-'
7492 write(ifhi,'(a,2i1)')'openhisto name ffom11',jaa,jex
7493 write(ifhi,'(a)')'htyp lba'
7494 write(ifhi,'(a)')'text 0.05 0.5 "+ ffom11a "'
7495 write(ifhi,'(a)')'array 2'
7498 z=ffom11a(u,-1.,jaa1,jaa2)
7499 write(ifhi,'(2e11.3)')u,z
7501 write(ifhi,'(a)') ' endarray'
7504 elseif(iii.eq.3)then
7506 write(ifhi,'(a)') '!----------------------------------'
7507 write(ifhi,'(a,3i1)') '! IB ',jaa,jex
7508 write(ifhi,'(a)') '!----------------------------------'
7510 !.......total integral
7521 z=zmin*(zmax/zmin)**(.5+tgss(ig1,i1)*(m1-1.5))
7524 if(xpmin.lt.xpmax)then
7527 xp=xpmin*(xpmax/xpmin)**(.5+tgss(ig2,i2)*(m2-1.5))
7529 r2=r2+wgss(ig2,i2)*ffsigiut(xp,xm,jaa,je1,je2)
7533 r2=r2*0.5*log(xpmax/xpmin)
7534 r1=r1+wgss(ig1,i1)*r2*z
7537 r1=r1*0.5*log(zmax/zmin)
7538 res= r1 * factk * .0390 /sigine*10
7541 xx2min = 0.01/engy !max(xpar1,0.01/engy)
7543 xx1min = 0.01/engy !max(xpar3,0.01/engy)
7545 nbins = 10 !nint(xpar5)
7547 write(ifhi,'(a,2i1)') 'openhisto xrange 0 1 name ffsig',jaa,jex
7548 write(ifhi,'(a)') 'yrange auto auto htyp lin xmod lin ymod log'
7549 write(ifhi,'(a)') 'txt "xaxis x+?IB! " '
7550 write(ifhi,'(a)') 'txt "yaxis dn?semi! / dx+?IB! "'
7551 write(ifhi,'(a,2i1,a)')'txt "title ffsig + MC (',jaa,jex,')"'
7552 write(ifhi,'(a)') 'array 2'
7553 del=(xx1max-xx1min)/nbins
7555 xx1=xx1min+(ii-0.5)*del
7560 xx2=xx2min*(xx2max/xx2min)**(.5+tgss(ig2,i2)*(m2-1.5))
7561 r2=r2+wgss(ig2,i2)*ffsigiut(xx1,xx2,jaa,je1,je2)*xx2
7564 sig=r2*0.5*log(xx2max/xx2min)
7565 sig = sig * factk * .0390 /sigine*10
7566 write(ifhi,'(2e12.4)')xx1,sig
7568 write(ifhi,'(a)') ' endarray'
7570 elseif(iii.eq.4)then
7572 write(ifhi,'(a)') '!----------------------------------'
7573 write(ifhi,'(a,3i1)') '! OB ',jaa,jex
7574 write(ifhi,'(a)') '!----------------------------------'
7584 pt=ptmin*(ptmax/ptmin)**(.5+tgss(ig,i)*(m-1.5))
7585 sig=ffsigi(pt**2,y2)
7586 sig =sig * factk * .0390 /sigine*10 * 2 ! 2 partons!
7587 sum=sum+wgss(ig,i)*sig*pt
7590 sum=sum*0.5*log(ptmax/ptmin)
7598 write(ifhi,'(a,i1)')'openhisto name jet',jj
7599 write(ifhi,'(a)')'xrange 0 20 xmod lin ymod log '
7600 write(ifhi,'(a)') 'txt "xaxis pt?OB! " '
7601 write(ifhi,'(a)') 'txt "yaxis dn?ptn! / dpt?OB! "'
7602 if(jj.eq.1)write(ifhi,'(a)')'htyp lro'
7603 if(jj.eq.2)write(ifhi,'(a)')'htyp lgo'
7604 if(jj.eq.3)write(ifhi,'(a)')'htyp lyo'
7605 write(ifhi,'(a,f7.2,a)') 'text 0.05 0.1 "1/f=',1./ff,'"'
7606 write(ifhi,'(a)')'array 2'
7607 delpt=(ptmax-ptmin)/nbins
7609 pt=ptmin+(i-0.5)*delpt
7612 sig=ffsigi(pt**2,y2) ! our stuff
7614 if(engy.ge.10.)sig=psjvrg1(pt**2,sx,y2) ! grv
7616 if(engy.ge.10.)sig=psjwo1(pt**2,sx,y2) !duke-owens
7618 sig =sig * factk * .0390 /sigine*10 * 2
7619 write(ifhi,'(2e12.4)')pt,sig
7621 write(ifhi,'(a)') ' endarray'
7622 if(jj.ne.1)write(ifhi,'(a)') 'closehisto'
7623 if(jj.ne.1)write(ifhi,'(a)') 'plot 0-'
7628 x=0.1+(min(3,iii)-2)*0.30
7629 y=0.2+(min(3,iii)-2)*0.55
7630 if(engy.gt.100.)then
7631 write(ifhi,'(a,2f5.2,a,f6.3,a)')'text',x,y,' " form ',sum,'"'
7633 write(ifhi,'(a,2f5.2,a,f6.5,a)')'text',x,y,' " form ',sum,'"'
7635 write(ifhi,'(a)') 'closehisto plot 0-'
7637 write(ifhi,'(a)') "!-----------------------------"
7638 write(ifhi,'(a)') "! MC "
7639 write(ifhi,'(a)') "!-----------------------------"
7642 * write(ifhi,'(a,i1,i1)')'openhisto name dndxPE',jaa,jex
7644 * write(ifhi,'(a,i1,i1)')'openhisto name dndxIB',jaa,jex
7646 * write(ifhi,'(a,i1,i1)')'openhisto name dndptOB',jaa,jex
7647 write(ifhi,'(a)') 'htyp prs'
7648 write(ifhi,'(a)') 'array 2'
7651 if(iii.eq.4)imax=nbipt
7654 if(iii.eq.4)u=ptob(i)
7659 if(iii.eq.2)z=z+wxe(kkk,jaai,jexi,i,j)
7660 if(iii.eq.3)z=z+wxb(kkk,jaai,jexi,i,j)
7661 if(iii.eq.4)z=z+wptob(jaai,jexi,i)
7666 if(iii.eq.4)del=(pto-ptu)/nbipt
7668 write(ifhi,'(2e11.3)')u,z
7671 write(ifhi,'(a)') ' endarray'
7672 x=0.1+(min(3,iii)-2)*0.30
7673 y=0.1+(min(3,iii)-2)*0.55
7675 write(ifhi,'(a,2f5.2,a,f6.3,a)')'text',x,y,' " simu ',sum,'"'
7677 write(ifhi,'(a,2f5.2,a,f6.5,a)')'text',x,y,' " simu ',sum,'"'
7679 write(ifhi,'(a)') 'closehisto'
7686 c-----------------------------------------------------------------------
7687 subroutine xEmsSe(iii,xmc,ptmc,ih,iqq)
7688 c-----------------------------------------------------------------------
7689 c iqq = 1 : String End mass and rapidity
7690 c iqq = 2 : String mass and rapidity
7691 c-----------------------------------------------------------------------
7696 common/cxpar/nx(2),x(nbix),wxmc(nbix,2),xmn,xmx,xu,xo
7698 common/cypar/ny(2),y(nbiy),wymc(nbiy,2),ymin,ymax,dy,yu,yo
7708 x(i)=xu*(xo/xu)**((i-0.5)/nbix)
7716 y(i)=yu+dy/2.+(i-1)*dy
7720 elseif(iii.eq.1)then
7723 if(ptmc.eq.0.)return
7724 if(iqq.eq.1)ymc=0.5*alog(xmc*s/ptmc)*ih
7725 if(iqq.eq.2)ymc=0.5*alog(xmc/ptmc)
7726 i=1+int(alog(xmc/xu)/alog(xo/xu)*nbix)
7729 wxmc(i,iqq)=wxmc(i,iqq)+1
7733 i=int((ymc-yu)/dy)+1
7736 wymc(i,iqq)=wymc(i,iqq)+1
7739 elseif(iii.eq.2)then
7741 write(ifhi,'(a)') '!--------------------------------'
7742 write(ifhi,'(a)') '! string end x distr '
7743 write(ifhi,'(a)') '!--------------------------------'
7744 write(ifhi,'(a)') 'openhisto'
7745 write(ifhi,'(a)') 'htyp lin'
7746 write(ifhi,'(a)') 'xmod log ymod log'
7747 write(ifhi,'(a,2e11.3)')'xrange',xu,xo
7748 if(iqq.eq.1)write(ifhi,'(a)') 'text 0 0 "xaxis string end x"'
7749 if(iqq.eq.2)write(ifhi,'(a)') 'text 0 0 "xaxis string x"'
7750 write(ifhi,'(a)') 'text 0 0 "yaxis P(x)"'
7751 write(ifhi,'(a)') 'array 2'
7753 dx=xu*(xo/xu)**(1.*i/nbix)*(1.-(xo/xu)**(-1./nbix))
7755 * write(ifhi,'(2e11.3)')x(i),wxmc(i,iqq)/dx/nx(iqq)
7757 write(ifhi,'(a)') ' endarray'
7758 write(ifhi,'(a)') 'closehisto plot 0'
7759 write(ifhi,'(a)') 'openhisto'
7760 write(ifhi,'(a)') 'htyp lin'
7761 write(ifhi,'(a)') 'xmod lin ymod lin'
7762 write(ifhi,'(a,2e11.3)')'xrange',yu,yo
7763 if(iqq.eq.1)write(ifhi,'(a)') 'text 0 0 "xaxis string end y"'
7764 if(iqq.eq.2)write(ifhi,'(a)') 'text 0 0 "xaxis string y"'
7765 write(ifhi,'(a)') 'text 0 0 "yaxis P(y)"'
7766 write(ifhi,'(a)') 'array 2'
7769 * write(ifhi,'(2e11.3)')y(i),wymc(i,iqq)/dy/ny(iqq)
7771 write(ifhi,'(a)') ' endarray'
7772 write(ifhi,'(a)') 'closehisto plot 0'
7778 c-----------------------------------------------------------------------
7779 subroutine xEmsDr(iii,xpmc,xmmc,ie)
7780 c-----------------------------------------------------------------------
7784 parameter(nbix=50,nie=4)
7785 common/cxpardr/nxp(nie),nxm(nie),x(nbix),wxpmc(nbix,nie)
7786 & ,wxmmc(nbix,nie),xmn,xmx,xu,xo,wxmc(nbix,nie),nx(nie)
7788 common/cypardr/ny(nie),y(nbiy),wymc(nbiy,nie),ymin,ymax,dy,yu,yo
7802 x(i)=xu*(xo/xu)**((i-0.5)/nbix)
7816 y(i)=yu+dy/2.+(i-1)*dy
7822 elseif(iii.eq.1)then
7824 if(ie.lt.1.or.ie.gt.nie)return
7826 if(xpmc.lt.xu)return
7827 i=1+int(alog(xpmc/xu)/alog(xo/xu)*nbix)
7830 wxpmc(i,ie)=wxpmc(i,ie)+1
7832 if(xmmc.lt.xu)return
7833 i=1+int(alog(xmmc/xu)/alog(xo/xu)*nbix)
7836 wxmmc(i,ie)=wxmmc(i,ie)+1
7840 ymc=0.5*alog(xpmc/xmmc)
7845 i=int((ymc-yu)/dy)+1
7848 wymc(i,ie)=wymc(i,ie)+1
7853 i=1+int(alog(xmc/xu)/alog(xo/xu)*nbix)
7856 wxmc(i,ie)=wxmc(i,ie)+1
7859 elseif(iii.eq.2)then
7863 if(ii.eq.1)write(ifhi,'(a)')'!----- projectile droplet ----'
7864 if(ii.eq.2)write(ifhi,'(a)')'!----- target droplet ----'
7865 if(ii.eq.3)write(ifhi,'(a)')'!----- projectile string end ----'
7866 if(ii.eq.4)write(ifhi,'(a)')'!----- target string end ----'
7867 write(ifhi,'(a)') '!--------------------------------'
7868 write(ifhi,'(a)') '! droplet/string x+ distr '
7869 write(ifhi,'(a)') '!--------------------------------'
7870 write(ifhi,'(a)') 'openhisto'
7871 write(ifhi,'(a)') 'htyp lru'
7872 write(ifhi,'(a)') 'xmod log ymod log'
7873 write(ifhi,'(a,2e11.3)')'xrange',xu,xo
7874 if(ii.eq.1.or.ii.eq.2)
7875 * write(ifhi,'(a)') 'text 0 0 "xaxis droplet x+"'
7876 if(ii.eq.3.or.ii.eq.4)
7877 * write(ifhi,'(a)') 'text 0 0 "xaxis string end x+"'
7878 write(ifhi,'(a)') 'text 0 0 "yaxis P(x)"'
7879 write(ifhi,'(a)') 'array 2'
7881 dx=xu*(xo/xu)**(1.*i/nbix)*(1.-(xo/xu)**(-1./nbix))
7883 * write(ifhi,'(2e11.3)')x(i),wxpmc(i,ii)/dx/nxp(ii)
7885 write(ifhi,'(a)') ' endarray'
7886 write(ifhi,'(a)') 'closehisto plot 0-'
7887 write(ifhi,'(a)') '!--------------------------------'
7888 write(ifhi,'(a)') '! droplet/string x- distr '
7889 write(ifhi,'(a)') '!--------------------------------'
7890 write(ifhi,'(a)') 'openhisto'
7891 write(ifhi,'(a)') 'htyp lba'
7892 write(ifhi,'(a)') 'xmod log ymod log'
7893 write(ifhi,'(a,2e11.3)')'xrange',xu,xo
7894 if(ii.eq.1.or.ii.eq.2)
7895 * write(ifhi,'(a)') 'text 0 0 "xaxis droplet x-"'
7896 if(ii.eq.3.or.ii.eq.4)
7897 * write(ifhi,'(a)') 'text 0 0 "xaxis string end x-"'
7898 write(ifhi,'(a)') 'text 0 0 "yaxis P(x)"'
7899 write(ifhi,'(a)') 'array 2'
7901 dx=xu*(xo/xu)**(1.*i/nbix)*(1.-(xo/xu)**(-1./nbix))
7903 * write(ifhi,'(2e11.3)')x(i),wxmmc(i,ii)/dx/nxm(ii)
7905 write(ifhi,'(a)') ' endarray'
7906 write(ifhi,'(a)') 'closehisto plot 0'
7907 write(ifhi,'(a)') '!--------------------------------'
7908 write(ifhi,'(a)') '! droplet/string y distr '
7909 write(ifhi,'(a)') '!--------------------------------'
7910 write(ifhi,'(a)') 'openhisto'
7911 write(ifhi,'(a)') 'htyp lin'
7912 write(ifhi,'(a)') 'xmod lin ymod lin'
7913 write(ifhi,'(a,2e11.3)')'xrange',yu,yo
7914 if(ii.eq.1.or.ii.eq.2)
7915 * write(ifhi,'(a)') 'text 0 0 "xaxis droplet y"'
7916 if(ii.eq.3.or.ii.eq.4)
7917 * write(ifhi,'(a)') 'text 0 0 "xaxis string end y"'
7918 write(ifhi,'(a)') 'text 0 0 "yaxis P(y)"'
7919 write(ifhi,'(a)') 'array 2'
7922 * write(ifhi,'(2e11.3)')y(i),wymc(i,ii)/dy/ny(ii)
7924 write(ifhi,'(a)') ' endarray'
7925 write(ifhi,'(a)') 'closehisto plot 0'
7929 write(ifhi,'(a)') '!--------------------------------'
7930 write(ifhi,'(a)') '! droplet/string mass distr '
7931 write(ifhi,'(a)') '!--------------------------------'
7935 if(ii.eq.2.or.ii.eq.4)write(ifhi,'(a)') 'closehisto plot 0-'
7936 if(ii.eq.3)write(ifhi,'(a)') 'closehisto plot 0'
7937 write(ifhi,'(a)') 'openhisto'
7938 if(ii.eq.1.or.ii.eq.3)write(ifhi,'(a)') 'htyp lru'
7939 if(ii.eq.2.or.ii.eq.4)write(ifhi,'(a)') 'htyp lba'
7940 write(ifhi,'(a)') 'xmod log ymod log'
7941 write(ifhi,'(a,2e11.3)')'xrange',sqrt(xu*s),sqrt(s*xo)
7942 if(ii.eq.1.or.ii.eq.2)
7943 * write(ifhi,'(a)') 'text 0 0 "xaxis droplet mass (GeV)"'
7944 if(ii.eq.4.or.ii.eq.3)
7945 * write(ifhi,'(a)') 'text 0 0 "xaxis string end mass (GeV)"'
7946 write(ifhi,'(a)') 'text 0 0 "yaxis P(x)"'
7947 write(ifhi,'(a)') 'array 2'
7949 dx=xu*(xo/xu)**(1.*i/nbix)*(1.-(xo/xu)**(-1./nbix))
7951 * write(ifhi,'(2e11.3)')sqrt(x(i)*s),wxmc(i,ii)/dx/nx(ii)
7953 write(ifhi,'(a)') ' endarray'
7955 write(ifhi,'(a)') 'closehisto plot 0'
7962 cc--------------------------------------------------------------------------
7963 c subroutine xtype(k,n,i1,i2,text)
7964 cc--------------------------------------------------------------------------
7966 c include 'epos.inc'
7967 c include 'epos.incems'
7968 c parameter(itext=40)
7973 c if(text(i:i).eq.'&')imax=i
7981 c write(ifch,*)('-',ll=1,27)
7982 c write(ifch,*)' '//text(1:imax-1)
7983 c write(ifch,*)('-',ll=1,27)
7988 c write(ifch,*)'k:',k,' n:',n,' ip:',ip,' it:',it
7989 c write(ifch,*)'bk:',bk(k)
7990 c if(n.ne.0)write(ifch,*)'idpr:',idpr(n,k)
7991 c write(ifch,*)'iep:',iep(ip),' iet:',iet(it)
7992 c write(ifch,*)'idp:',idp(ip),' idt:',idt(it)
7997 c------------------------------------------------------------------------
7998 subroutine XPrint(text)
7999 c------------------------------------------------------------------------
8001 include 'epos.incems'
8002 double precision xpptot,xmptot,xpttot,xmttot
8007 if(text(i:i).eq.'&')imax=i
8009 write(ifch,'(1x,a)')text(1:imax-1)
8011 write(ifch,'(a)')' npr0: npr1: nprmx: Pomeron id lattice:'
8013 write(ifch,'(1x,i6,1x,i2,6x,i2,6x,i2,7x,$)')
8014 * k,npr(0,k),npr(1,k),nprmx(k)
8016 write(ifch,'(i2,$)')idpr(n,k)
8025 write(ifch,'(a)')' Pomeron xy lattice:'
8028 xpptot=xpptot+xppr(n,k)
8029 xmttot=xmttot+xmpr(n,k)
8030 write(ifch,'(i6,1x,i2,1x,d10.3,1x,d10.3,3x,$)')
8031 * k,n,xpr(n,k),ypr(n,k)
8036 write(ifch,'(a)')' projectile remnants x+,x-,px,py,x,iep:'
8038 xpptot=xpptot+xpp(ip)
8039 xmptot=xmptot+xmp(ip)
8040 write(ifch,'(i3,2x,5d12.3,i3)')ip,xpp(ip),xmp(ip),xxp(ip),xyp(ip)
8044 write(ifch,'(a)')' target remnants x-,x+,px,py,x,iet:'
8046 xpttot=xpttot+xpt(it)
8047 xmttot=xmttot+xmt(it)
8048 write(ifch,'(i3,2x,5d12.3,i3)')it,xmt(it),xpt(it),xxt(it),xyt(it)
8052 write(ifch,*)' remnant balance x+,x-:'
8053 &,(xpptot+xpttot)/dble(maproj)
8054 &,(xmptot+xmttot)/dble(matarg)
8058 c-------------------------------------------------------------------------
8060 c-------------------------------------------------------------------------
8062 double precision fom,x
8063 write(ifhi,'(a)') '!##################################'
8064 write(ifhi,'(a,i3)') '! fom '
8065 write(ifhi,'(a)') '!##################################'
8069 xi=0.01+0.16*float(i-1)
8070 write(ifhi,'(a,i1)') 'openhisto name fom',i
8071 write(ifhi,'(a)') 'htyp lin xmod lin ymod log'
8072 write(ifhi,'(a)') 'xrange 0 1'
8073 write(ifhi,'(a)') 'yrange 0.1 1000 '
8074 write(ifhi,'(a)') 'text 0 0 "xaxis x "'
8075 write(ifhi,'(a)') 'text 0 0 "yaxis fom"'
8077 & write(ifhi,'(a,f4.2,a,f4.1,a)')'text ',xi,' 0.9 "',z,'"'
8079 & write(ifhi,'(a,f4.2,a,f4.0,a)')'text ',xi,' 0.9 "',z,'"'
8080 write(ifhi,'(a)') 'array 2'
8083 write(ifhi,'(2e11.3)')x,fom(z,x,b)
8085 write(ifhi,'(a)') ' endarray'
8086 write(ifhi,'(a)') ' closehisto '
8087 if(i.lt.6)write(ifhi,'(a)') 'plot 0-'
8088 if(i.eq.6)write(ifhi,'(a)') 'plot 0'