1 c-----------------------------------------------------------------------
2 subroutine utresc(iret)
3 c-----------------------------------------------------------------------
4 c if irescl=1 rescaling is done, otherwise the purpose of going through
5 c this routine is to not change the seed in case of no rescaling
6 c-----------------------------------------------------------------------
9 double precision p1,esoll,ppp,seedp
10 dimension p1(5),p0(5,mamx+mamx),pini(mxptl)
11 logical force,nolead(mxptl),lspec(mxptl)
14 call utpri('utresc',ish,ishini,4)
16 errlim=0.001 !max(0.001,1./engy)
20 nptlpt=iabs(maproj)+iabs(matarg)
21 call ranfgt(seedp) !not to change the seed ...
22 if(nptl.le.nptlpt) goto 9999
40 c if(mod(istptl(i),10).eq.1)then
42 p1(j)=p1(j)+dble(pptl(j,i))
47 if(mod(istptl(i),10).eq.0)then
49 p1(j)=p1(j)+dble(pptl(j,i))
52 if((ityptl(i).eq.45.and.maproj.ge.100)
53 & .or.(ityptl(i).eq.55.and.matarg.ge.100))lspec(i)=.true.
54 if((abs(pptl(3,i)/pnullx).le.0.9
55 & .and.abs(pptl(3,i)).gt.pptl(5,i)).or.lspec(i))then
57 c write(ifch,*)'nolead',i
60 c write(ifch,*)'lead',i
64 ppp=(p1(4)+p1(3))*(p1(4)-p1(3))-p1(2)*p1(2)-p1(1)*p1(1)
69 write(ifch,*)'p1',p1(1),p1(2),p1(3),p1(4),ppp
70 if(ish.ge.2)write (ifch,*) 'Problem in utresc (0): redo event'
71 write (ifmt,*) 'Problem in utresc (0): redo event'
72 c call utstop('utresc&')
76 if(ish.ge.4) write (ifch,'(a,5g13.6)') 'boost-vector: ',p1
82 if(mod(istptl(i),10).le.1)then
83 call utlob4(1,p1(1),p1(2),p1(3),p1(4),p1(5)
84 $ ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i))
86 if(i.gt.nptlpt.and.nolead(i))pmax=max(pmax,abs(pptl(3,i)))
90 call alistf('list after boost&')
93 if(ish.ge.5)write(ifch,'(a)')'--------rescale momenta----------'
95 c rescale momenta in rest frame
96 c -----------------------------
108 nfirst=int(rangen()*float(npart)) !start list randomly
111 if(j.gt.nptl)j=imin+i+nfirst-npart
112 if(mod(istptl(j),10).eq.0)then
115 c if(j.gt.nptlpt)then
116 c if(abs(pptl(3,j))/pnullx.lt.0.9)then !not spectator or diffraction
117 if(scal.eq.1..and.diff0.eq.0.)then
123 & ityptl(j)/10.eq.4.or.ityptl(j)/10.eq.5
124 & ))then !try just remnant first
126 diff=sign(min(0.3*abs(pini(j)),
127 & rangen()*abs(difft)),difft)
128 pptl3new=scal*(pptl(3,j)-diff)
129 c write(ifch,*)'par',j,pptl3new,pptl(3,j),diff,difft
131 if(abs(pptl3new).lt.pmax)then
132 if(abs(pptl3new-pini(j)).lt.ferr*abs(pini(j))
133 & .or.(lspec(j).and.abs(pptl3new).lt.abs(0.8*pini(j))))then !particle should not be too fast or too modified
134 c write(ifch,*)'used'
136 pptl(3,j)=scal*(pptl(3,j)-diff)
137 pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2
138 * +pptl(3,j)**2+pptl(5,j)**2)
144 c sum over all particles
151 scal=real(esoll)/(sum-diff)
152 if(ish.ge.6)write(ifch,*)
153 $ 'ipass,scal,diff/pnullx,e,esoll,pz,pmax,ndif,f:'
154 $ ,ipass,scal,diff/pnullx,sum,esoll,sum3,pnullx,ndif,force
155 if(abs(scal-1.).le.errlim.and.abs(diff/pnullx).lt.errlim)
157 if(ndif.gt.0.and.(force.or.ipass.lt.150))then
160 elseif(abs(scal-1.).le.1e-2.and.abs(diff/pnullx).lt.5e-2)then
164 $ write(ifmt,*)'Warning in utresc: no more allowed particle'
173 call utmsg('utrescl')
174 write(ifch,*)'***** scal=',scal,diff
177 if(abs(scal)+abs(diff/pnullx).gt.2.5)then
178 write(ifmt,*)'Warning in utresc !'
179 write(ifch,*)'redo EPOS event ...'
193 if(mod(istptl(i),10).le.1)then
194 call utlob4(-1,p1(1),p1(2),p1(3),p1(4),p1(5)
195 $ ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i))
200 if(ish.ge.4)call alist('list after rescaling&',1,nptl)
204 scalmean=scalmean+scal
205 write(ifch,*)' average rescaling factor: ',scalmean
208 call ranfst(seedp) ! ... after this subroutine
209 call utprix('utresc',ish,ishini,4)
213 c-----------------------------------------------------------------------
214 subroutine utghost(iret)
215 c-----------------------------------------------------------------------
216 c if irescl=1 make particle on-shell if not
217 c-----------------------------------------------------------------------
219 include 'epos.incems'
220 double precision seedp
221 call utpri('ughost',ish,ishini,4)
224 nptlpt=iabs(maproj)+iabs(matarg)
225 call ranfgt(seedp) !not to change the seed ...
226 if(nptl.le.nptlpt) goto 9999
228 if(ish.ge.5)write(ifch,'(a)')'---------mark ghosts---------'
233 if(mod(istptl(j),10).eq.0.and.pptl(4,j).gt.0.d0)then
235 call idmass(idptl(j),amass)
236 if(abs(pptl(5,j)-amass).gt.0.01*amass)then
237 if(ish.ge.5)write(ifch,*)'wrong particle mass',j,idptl(j)
240 call idres(idptl(j),amass,idr,iadj,0)
245 call idmass(idptl(j),amass)
248 call idtau(idptl(j),pptl(4,j),pptl(5,j),taugm)
249 tivptl(2,j)=tivptl(1,j)+taugm*(-alog(rangen()))
251 if(abs((pptl(4,j)+pptl(3,j))*(pptl(4,j)-pptl(3,j))
252 $ -pptl(2,j)**2-pptl(1,j)**2-pptl(5,j)**2).gt.0.3
253 $ .and.abs(1.-abs(pptl(3,j))/pptl(4,j)).gt.0.01)then
254 !print*,'ghost',ityptl(j),idptl(j)
255 if(ish.ge.1)write(ifmt,*)'ghost:',j,idptl(j),ityptl(j)
257 write(ifch,'(a,$)')'ghost:'
260 ityptl(j)=100+ityptl(j)/10
262 elseif(mod(istptl(j),10).eq.0)then
264 write(ifmt,*)'Lost particle (E=0)'
265 write(ifch,*)'Lost particle (E=0) :'
266 call alistc("utghost&",j,j)
268 istptl(j)=istptl(j)+2
272 if(ish.ge.5)write(ifch,'(a)')'---------treat ghosts---------'
286 if(mod(istptl(j),10).eq.0)then
287 if(ityptl(j).ge.101.and.ityptl(j).le.105)then
291 if(pptl(4,j).gt.0.)efif=efif+pptl(4,j)
295 if(pptl(4,j).gt.0.)then
297 amt=1.-(pptl(5,j)*Einv)**2+(pptl(1,j)*Einv)**2
298 $ +(pptl(2,j)*Einv)**2
303 pptl(3,j)=sign(pptl(4,j),pptl(3,j))*sqrt(amt)
305 y=(rangen()+rangen()+rangen()+rangen()-2.)/2.*yhaha
306 y=sign(abs(y),pptl(3,j))
308 $ =sqrt(pptl(5,j)**2+pptl(1,j)**2
309 $ +pptl(2,j)**2)*sinh(y)
311 $ =sqrt(pptl(5,j)**2+pptl(1,j)**2
312 $ +pptl(2,j)**2)*cosh(y)
319 pptl(k,j)=pptl(k,j)*scal
321 pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2
328 $ write (ifch,*) 'nrevt,psum,esum,pfif,efif,nfif,scal'
329 $ ,nrevt,psum,esum,pfif,efif,nfif,scal
333 if ( ish.ge.5 ) write (ifch,*) 'tot',nfif,efif,pfif,esum,psum
336 if(nfif.gt.5.or.(esum.gt.0.05*engy.and.nfif.ne.1))then
339 if ( ityptl(j).ge.101 .and. ityptl(j).le.105 )then
340 if((psum-pfif)*(1.-scal).ge.0)
341 & pptl(3,j)=pptl(3,j)-(psum-pfif)/nfif
349 if ( ish.ge.5 ) write (ifch,*) 'scal',scal
350 if ( abs(scal-1.) .gt. 0.05 ) then
355 if(ish.ge.2)write (ifch,*) 'Problem in utghost : redo event'
356 if(ish.ge.1)write (ifmt,*) 'Problem in utghost : redo event'
362 if ( ityptl(j).ge.101 .and. ityptl(j).le.105 )then
363 pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2
369 if(ish.ge.5)write(ifch,'(a)')'---------Check Ghost list---------'
375 if(mod(istptl(j),10).eq.0)then
376 if(ityptl(j).le.105.and.ityptl(j).ge.101)then
377 write(ifch,'(a,$)')'ghost:'
386 call ranfst(seedp) ! ... after this subroutine
387 call utprix('ughost',ish,ishini,4)
391 c-----------------------------------------------------------------------
392 subroutine utrsph(iret)
393 c-----------------------------------------------------------------------
394 c if irescl=1 and ispherio=1 rescaling is done for particle used by
395 c spherio as initial condition.
396 c-----------------------------------------------------------------------
398 include 'epos.incems'
399 double precision p1,esoll
400 dimension p1(5),p0(5,mamx+mamx)
401 call utpri('utrsph',ish,ishini,4)
406 nptlpt=iabs(maproj)+iabs(matarg)
407 if(nptl.le.nptlpt) goto 9999
416 $ .and.(iorptl(i).ge.1.and.istptl(iorptl(i)).eq.41))
417 $ .or.istptl(i).eq.20.or.istptl(i).eq.21)then
419 p1(j)=p1(j)+dble(pptl(j,i))
428 p1(j)=p1(j)+dble(pptl(j,i))
431 p1(5)=dsqrt((p1(4)+p1(3))*(p1(4)-p1(3))-p1(2)**2.d0-p1(1)**2.d0)
433 if(ish.ge.4) write (ifch,'(a,5g13.6)') 'boost-vector',p1
439 $ .and.(iorptl(i).ge.1.and.istptl(iorptl(i)).eq.41))
440 $ .or.istptl(i).eq.20.or.istptl(i).eq.21
441 $ .or.(istptl(i).eq.0.and.i.le.nptlpt))then
442 call utlob4(1,p1(1),p1(2),p1(3),p1(4),p1(5)
443 $ ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i))
448 if(ish.ge.5)write(ifch,'(a)')'------------------'
450 c rescale momenta in rest frame
451 c -----------------------------
461 $ .and.(iorptl(j).ge.1.and.istptl(iorptl(j)).eq.41))
462 $ .or.istptl(j).eq.20.or.istptl(j).eq.21
463 $ .or.(istptl(j).eq.0.and.j.le.nptlpt))then
466 pptl(3,j)=scal*(pptl(3,j)-diff)
467 pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2
477 if(ish.ge.6)write(ifch,*)'ipass,scal,diff,e,esoll,pz,ndif:'
478 $ ,ipass,scal,diff,sum,esoll,sum3,ndif
479 if(abs(scal-1.).le.errlim.and.abs(diff).lt.10.*errlim) goto300
483 write(ifch,*)'***** scal=',scal,diff
494 $ .and.(iorptl(i).ge.1.and.istptl(iorptl(i)).eq.41))
495 $ .or.istptl(i).eq.20.or.istptl(i).eq.21
496 $ .or.(istptl(i).eq.0.and.i.le.nptlpt))then
497 call utlob4(-1,p1(1),p1(2),p1(3),p1(4),p1(5)
498 $ ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i))
507 9999 call utprix('utrsph',ish,ishini,4)
511 cc-----------------------------------------------------------------------
512 c double precision function dddlog(xxx)
513 cc-----------------------------------------------------------------------
514 c double precision xxx
516 c if(xxx.gt.0d0)dddlog=dlog(xxx)
519 ccc-----------------------------------------------------------------------
520 c subroutine randfl(jc,iqa0,iflav,ic,isame)
521 cc-----------------------------------------------------------------------
522 cc returns random flavour ic(2) (iqa0=1:quark,2:antiquark,11:diquark)
523 cc-----------------------------------------------------------------------
525 c real probab(nflav),probsu(nflav+1)
526 c integer jc(nflav,2),jc0(nflav,2),ic(2)
528 c write(ifch,*)('-',i=1,10)
529 c *,' entry sr randfl ',('-',i=1,30)
530 c write(ifch,*)'iqa0:',iqa0
542 c if(iqa1.eq.0)goto9998
546 c probab(i)=jc(i,iqa)-jc0(i,iqa)
547 c if(isame.eq.1)probab(i)=probab(i)*(jc(i,3-iqa)-jc0(i,3-iqa))
557 c probsu(i+1)=probsu(i)+probab(i)/su
558 c if(probsu(i+1)-probsu(i).lt.1e-5)probsu(i+1)=probsu(i)
560 c r=rangen()*probsu(nflav+1)
562 c if(probsu(i).le.r.and.r.lt.probsu(i+1))iflav=i
564 c jc0(iflav,iqa)=jc0(iflav,iqa)+1
565 c if(isame.eq.1)jc0(iflav,3-iqa)=jc0(iflav,3-iqa)+1
566 c call idenco(jc0,ic,ireten)
567 c if(ireten.eq.1)call utstop('randfl: idenco ret code = 1&')
569 c write(ifch,*)'probab:'
570 c write(ifch,*)probab
571 c write(ifch,*)'probsu:'
572 c write(ifch,*)probsu
573 c write(ifch,*)'ran#:',r,' flav:',iflav
577 c if(ish.ge.6)write(ifch,*)('-',i=1,30)
578 c *,' exit sr randfl ',('-',i=1,10)
583 cc-----------------------------------------------------------------------
584 c subroutine ranhvy(x,eps)
585 cc-----------------------------------------------------------------------
586 cc generates x for heavy particle fragmentation according to
588 cc d(x)=1/(x*(1-1/x-eps/(1-x))**2)
589 cc =d0(x)*d1(x)*d2(x)
590 cc d0(x)=(1-x)**2/((1-x)**2+eps)**2
592 cc d2(x)=(((1-x)**2+eps)/((1-x)**2+eps*x))**2
594 cc generates flat in x if eps>1.
595 cc-----------------------------------------------------------------------
598 c pow=alog((3.+eps)/eps)/aln4
599 c ymx=(eps*(3.*pow-1.)/(pow+1.))**(.5/pow)
601 c d0mx=(1-zmx)**2/((1.-zmx)**2+eps)**2*pow*ymx**(pow-1.)
602 c d2mx=2./(2.-sqrt(eps))
606 c d0mx=(1.-zmx)**2/((1.-zmx)**2+eps)**2
610 cc generate z according to (1-z)**2/((1-z)**2+eps*z)**2
615 c d0z=(1.-z)**2/((1.-z)**2+eps)**2*pow*y**(pow-1.)
616 c if(d0z.lt.rangen()*d0mx) goto1
618 cc check remaining factors
620 c d2=(((1.-z)**2+eps)/((1.-z)**2+eps*z))**2
621 c if(d1*d2.lt.rangen()*d2mx) goto1
628 c-----------------------------------------------------------------------
630 c-----------------------------------------------------------------------
631 c returns randomly +1 or -1
632 c-----------------------------------------------------------------------
634 if(rangen().gt.0.5)ransig=-1
638 cc-----------------------------------------------------------------------
639 c function ranxq(n,x,q,xmin)
640 cc-----------------------------------------------------------------------
641 cc returns random number according to x(i) q(i) with x>=xmin
642 cc-----------------------------------------------------------------------
646 c if(xmin.eq.0.)goto3
650 c if(x(i).lt.xmin)then
652 c elseif(x(i).gt.xmin)then
658 c if(i2-i1.gt.1)goto1
661 c if(q(imin).gt.q(n)*.9999)then
665 c qran=q(imin)+rangen()*(q(n)-q(imin))
666 c ranxq=utinvt(n,x,q,qran)
669 c if(ranxq.lt.xmin)then
670 c call utmsg('ranxq ')
671 c write(ifch,*)'***** ranxq=',ranxq,' < xmin=',xmin
672 c write(ifch,*)'q(imin) q q(n):',q(imin),qran,q(n)
673 c write(ifch,*)'x(imin) x x(n):',x(imin),ranxq,x(n)
681 cc ***** end r-routines
682 cc ***** beg s-routines
684 cc-----------------------------------------------------------------------
686 cc-----------------------------------------------------------------------
687 c sbet=utgam1(z)*utgam1(w)/utgam1(z+w)
691 cc-----------------------------------------------------------------------
692 c function smass(a,y,z)
693 cc-----------------------------------------------------------------------
694 cc returns droplet mass (in gev) (per droplet, not (!) per nucleon)
695 cc according to berger/jaffe mass formula, prc35(1987)213 eq.2.31,
696 cc see also c. dover, BNL-46322, intersections-meeting, tucson, 91.
697 cc a: massnr, y: hypercharge, z: charge,
698 cc-----------------------------------------------------------------------
699 c common/cmass/thet,epsi,as,ac,dy,dz,ym,cz,zm,sigma,rzero
701 c zmin=cz/(dz/a+zm/a**(1./3.))
702 c smass=epsi*a+as*a**(2./3.)+(ac/a**(1./3.)+dz/a/2.)*(z-zmin)**2
703 c *+dy/a/2.*(y-ymin)**2
707 cc-----------------------------------------------------------------------
708 c subroutine smassi(theta)
709 cc-----------------------------------------------------------------------
710 cc initialization for smass.
711 cc calculates parameters for berger/jaffe mass formula
712 cc (prc35(1987)213 eq.2.31, see also c. dover, BNL-46322).
713 cc theta: parameter that determines all parameters in mass formula.
714 cc-----------------------------------------------------------------------
715 c common/cmass/thet,epsi,as,ac,dy,dz,ym,cz,zm,sigma,rzero
725 c rzero=si/astr/( 2./3./pi*(1+co**3) )**(1./3.)
726 cctp060829 cs=astr/si
727 c cz=-astr/si*( ( .5*(1+co**3) )**(1./3.)-1 )
728 c sigma=6./8./pi*(astr/si)**3*(co**2/6.-si**2*(1-si)/3.-
729 c *1./3./pi*(pi/2.-theta-sin(2*theta)+si**3*alog((1+co)/si)))
731 c epsi=astr*((.5*(1+co**3))**(1./3.)+2)/si
732 c as=4*pi*sigma*rzero**2
734 c dz=astr/si*bet**(1./3.)*co**2*
735 c *(co**4*(1+bet**(2./3.))+(1+bet)**2)/
736 c *( (2*co**2+bet**(1./3.))*(co**4*(1+bet**(2./3.))+(1+bet)**2)-
737 c *(co**4+bet**(1./3.)*(1+bet))*((2*bet**(2./3.)-1)*co**2+1+bet) )
738 c dy=astr/6.*(1+co**3)**3/si*
739 c *( 1+(1+co)/(4*(1+co**3))**(2./3.) )/
740 c *(co**6+co+co*(.5*(1+co**3))**(4./3.))
742 c ym=(1-co**3)/(1+co**3)
747 cc-----------------------------------------------------------------------
749 cc-----------------------------------------------------------------------
751 cc-----------------------------------------------------------------------
753 c common/cmass/thet,epsi,as,ac,dy,dz,ym,cz,zm,sigma,rzero
754 c real eng(14),ymi(14),zmi(14)
756 c write(ifch,*)'parameters of mass formula:'
757 c write(ifch,*)'theta=',thet,' epsi=',epsi
758 c write(ifch,*)'as=',as,' ac=',ac
759 c write(ifch,*)'dy=',dy,' dz=',dz
760 c write(ifch,*)'ym=',ym
761 c write(ifch,*)'cz dz zm=',cz,dz,zm
762 c write(ifch,*)'sigma**1/3=',sigma**(1./3.),' rzero=',rzero
763 c write(ifch,*)'mass:'
764 c write(ifch,5000)(j,j=1,14)
765 c5000 format(5x,'a:',14i5)
769 c4 zmi(j)=cz/(dz/a+zm/a**(1./3.))
770 c write(ifch,5002)(ymi(j),j=1,14)
771 c5002 format(1x,'ymin: ',14f5.2)
772 c write(ifch,5003)(zmi(j),j=1,14)
773 c5003 format(1x,'zmin: ',14f5.2)
780 c3 eng(j)=smass(a,y,z)/a
781 c write(ifch,5001)ns,(eng(j),j=1,14)
782 c5001 format(1x,'s=',i2,2x,14f5.2)
784 c write(ifch,*)'mass-mass(free):'
785 c write(ifch,5000)(j,j=1,14)
792 c call smassu(a,y,z,ku,kd,ks,kc)
793 c6 eng(j)=(smass(a,y,z)-utamnu(ku,kd,ks,kc,0,0,3))/a
794 c write(ifch,5001)ns,(eng(j),j=1,14)
800 cc-----------------------------------------------------------------------
801 c subroutine smasst(kux,kdx,ksx,kcx,a,y,z)
802 cc-----------------------------------------------------------------------
803 cc input: kux,kdx,ksx,kcx = net quark numbers (for u,d,s,c quarks).
804 cc output: massnr a, hypercharge y and charge z.
805 cc-----------------------------------------------------------------------
807 c if(kux+kdx+ksx+kcx.lt.0.)sg=-1
813 c if(mod(k,3).ne.0)stop'noninteger baryon number'
817 c if(mod(nz,3).ne.0)stop'noninteger charge'
822 cc-----------------------------------------------------------------------
823 c subroutine smassu(ax,yx,zx,ku,kd,ks,kc)
824 cc-----------------------------------------------------------------------
825 cc input: massnr ax, hypercharge yx and charge zx.
826 cc output: ku,kd,ks,kc = net quark numbers (for u,d,s,c quarks).
827 cc-----------------------------------------------------------------------
840 cc-----------------------------------------------------------------------
841 c function spoc(a,b,c,d,x)
842 cc-----------------------------------------------------------------------
843 cc power fctn with cutoff
844 cc-----------------------------------------------------------------------
846 c if(a.eq.0..and.b.eq.0.)return
849 c spoc=amin1(spoc,spoc0)
850 c spoc=amax1(0.,spoc)
854 c-----------------------------------------------------------------------
856 c-----------------------------------------------------------------------
857 c returns acos(x) for -1 <= x <= 1 , acos(+-1) else
858 c-----------------------------------------------------------------------
864 write(ifch,*)'***** argum = ',argum,' set -1'
871 write(ifch,*)'***** argum = ',argum,' set 1'
880 c----------------------------------------------------------------------
881 function utamnu(keux,kedx,kesx,kecx,kebx,ketx,modus)
882 c----------------------------------------------------------------------
883 c returns min mass of droplet with given u,d,s,c content
884 c keux: net u quark number
885 c kedx: net d quark number
886 c kesx: net s quark number
887 c kecx: net c quark number
888 c kebx: net b quark number
889 c ketx: net t quark number
890 c modus: 4=two lowest multiplets; 5=lowest multiplet
891 c----------------------------------------------------------------------
892 common/files/ifop,ifmt,ifch,ifcx,ifhi,ifdt,ifcp,ifdr
893 common/csjcga/amnull,asuha(7)
894 common/drop4/asuhax(7),asuhay(7)
896 if(modus.lt.4.or.modus.gt.5)stop'UTAMNU: not supported'
897 1 format(' flavours:',6i5 )
898 100 format(' flavours+mass:',6i5,f8.2 )
899 c write(ifch,1)keux,kedx,kesx,kecx,kebx,ketx c write(ifmt,*)'wrong mass in gethad ',damss
904 if(modus.eq.4)asuha(i)=asuhax(i) !two lowest multiplets
905 if(modus.eq.5)asuha(i)=asuhay(i) !lowest multiplet
908 ke=iabs(keux+kedx+kesx+kecx+kebx+ketx)
910 if(keux+kedx+kesx+kecx+kebx+ketx.ge.0)then
926 c write(ifch,*)keu,ked,kes,kec,keb,ket
928 c removing top mesons to remove t quarks or antiquarks
933 if(ii*keu.le.ii*ked)then
938 amnull=amnull+200. ! ???????
942 c removing bottom mesons to remove b quarks or antiquarks
947 if(ii*keu.le.ii*ked)then
952 amnull=amnull+5.28 ! (B-meson)
956 c removing charm mesons to remove c quarks or antiquarks
961 if(keu*ii.le.ked*ii)then
966 amnull=amnull+1.87 ! (D-meson)
970 c write(ifch,100)keu,ked,kes,kec,keb,ket,amnull
972 c removing mesons to remove s antiquarks
975 amnull=amnull+asuha(6)
985 c removing mesons to remove d antiquarks
989 amnull=amnull+asuha(5)
992 amnull=amnull+asuha(6)
999 c removing mesons to remove u antiquarks
1003 amnull=amnull+asuha(5)
1006 amnull=amnull+asuha(6)
1013 c write(ifch,100)keu,ked,kes,kec,keb,ket,amnull
1014 c print*,keu,ked,kes,kec,keb,ket,amnull
1016 if(keu+ked+kes+kec+keb+ket.ne.ke)
1017 *call utstop('utamnu: sum_kei /= ke&')
1022 c removing strange baryons
1026 if((4-i)*kes.gt.(i-1)*keq)then
1027 amnux=amnux+asuha(1+i)
1030 if(kes.lt.0)call utstop('utamnu: negative kes&')
1031 if(keq.lt.0)call utstop('utamnu: negative keq&')
1045 if(keu+ked.ne.keq)call utstop('utamnu: keu+ked /= keq&')
1046 c write(ifch,100)keu,ked,kes,kec,keb,ket,amnull+amnux
1047 c print*,keu,ked,kes,kec,keb,ket,amnull+amnux
1049 c removing nonstrange baryons
1051 if(keu.gt.2*ked)then
1052 amnux=amnux+asuha(7)
1054 if(keu.lt.0)call utstop('utamnu: negative keu&')
1057 if(ked.gt.2*keu)then
1058 amnux=amnux+asuha(7)
1060 if(ked.lt.0)call utstop('utamnu: negative ked&')
1065 c write(ifch,100)keu,ked,kes,kec,keb,ket,amnull+amnux
1066 c print*,keu,ked,kes,kec,keb,ket,amnull+amnux
1068 if(mod(keq,3).ne.0)call utstop('utamnu: mod(keq,3) /= 0&')
1069 amnux=amnux+asuha(1)*keq/3
1071 c write(ifch,100)keu,ked,kes,kec,keb,ket,amnull+amnux
1072 c print*,keu,ked,kes,kec,keb,ket,amnull+amnux
1076 if(amnull.eq.0)amnull=asuha(5)
1082 c-----------------------------------------------------------------------
1083 function utamnx(jcp,jcm)
1084 c-----------------------------------------------------------------------
1085 c returns minimum mass for the decay of jcp---jcm (by calling utamnu).
1086 c-----------------------------------------------------------------------
1088 integer jcp(nflav,2),jcm(nflav,2)
1092 if(jcp(i,j).ne.0)goto1
1095 keu=jcm(1,1)-jcm(1,2)
1096 ked=jcm(2,1)-jcm(2,2)
1097 kes=jcm(3,1)-jcm(3,2)
1098 kec=jcm(4,1)-jcm(4,2)
1099 keb=jcm(5,1)-jcm(5,2)
1100 ket=jcm(6,1)-jcm(6,2)
1101 utamnx=utamnu(keu,ked,kes,kec,keb,ket,5)
1107 if(jcm(i,j).ne.0)goto2
1110 keu=jcp(1,1)-jcp(1,2)
1111 ked=jcp(2,1)-jcp(2,2)
1112 kes=jcp(3,1)-jcp(3,2)
1113 kec=jcp(4,1)-jcp(4,2)
1114 keb=jcp(5,1)-jcp(5,2)
1115 ket=jcp(6,1)-jcp(6,2)
1116 utamnx=utamnu(keu,ked,kes,kec,keb,ket,5)
1120 keu=jcp(1,1)-jcp(1,2)
1121 ked=jcp(2,1)-jcp(2,2)
1122 kes=jcp(3,1)-jcp(3,2)
1123 kec=jcp(4,1)-jcp(4,2)
1124 keb=jcp(5,1)-jcp(5,2)
1125 ket=jcp(6,1)-jcp(6,2)
1126 ke=keu+ked+kes+kec+keb+ket
1127 if(mod(ke+1,3).eq.0)then
1129 amms1=utamnu(keu,ked,kes,kec,keb,ket,5)
1132 amms2=utamnu(keu,ked,kes,kec,keb,ket,5)
1133 elseif(mod(ke-1,3).eq.0)then
1135 amms1=utamnu(keu,ked,kes,kec,keb,ket,5)
1138 amms2=utamnu(keu,ked,kes,kec,keb,ket,5)
1140 call utstop('utamnx: no singlet possible (1)&')
1142 keu=jcm(1,1)-jcm(1,2)
1143 ked=jcm(2,1)-jcm(2,2)
1144 kes=jcm(3,1)-jcm(3,2)
1145 kec=jcm(4,1)-jcm(4,2)
1146 keb=jcm(5,1)-jcm(5,2)
1147 ket=jcm(6,1)-jcm(6,2)
1148 ke=keu+ked+kes+kec+keb+ket
1149 if(mod(ke+1,3).eq.0)then
1151 amms3=utamnu(keu,ked,kes,kec,keb,ket,5)
1154 amms4=utamnu(keu,ked,kes,kec,keb,ket,5)
1155 elseif(mod(ke-1,3).eq.0)then
1157 amms3=utamnu(keu,ked,kes,kec,keb,ket,5)
1160 amms4=utamnu(keu,ked,kes,kec,keb,ket,5)
1162 call utstop('utamnx: no singlet possible (2)&')
1164 utamnx=min(amms1+amms3,amms2+amms4)
1165 c print *,amms1,amms3,amms2,amms4,jcp,jcm
1171 cc-----------------------------------------------------------------------
1172 c function utamny(jcp,jcm)
1173 cc-----------------------------------------------------------------------
1174 cc returns minimum mass of jcp+jcm (by calling utamnu).
1175 cc-----------------------------------------------------------------------
1176 c parameter (nflav=6)
1177 c integer jcp(nflav,2),jcm(nflav,2),jc(nflav,2)
1179 c jc(nf,1)=jcp(nf,1)+jcm(nf,1)
1180 c7 jc(nf,2)=jcp(nf,2)+jcm(nf,2)
1181 c keu=jc(1,1)-jc(1,2)
1182 c ked=jc(2,1)-jc(2,2)
1183 c kes=jc(3,1)-jc(3,2)
1184 c kec=jc(4,1)-jc(4,2)
1185 c keb=jc(5,1)-jc(5,2)
1186 c ket=jc(6,1)-jc(6,2)
1187 c utamny=utamnu(keu,ked,kes,kec,keb,ket,5)
1191 c-----------------------------------------------------------------------
1192 function utamnz(jc,modus)
1193 c-----------------------------------------------------------------------
1194 c returns minimum mass of jc (by calling utamnu).
1195 c-----------------------------------------------------------------------
1204 utamnz=utamnu(keu,ked,kes,kec,keb,ket,modus)
1208 c-----------------------------------------------------------------------
1209 subroutine utar(i1,i2,i3,x0,x1,x2,x3,xx)
1210 c-----------------------------------------------------------------------
1211 c returns the array xx with xx(1)=x0 <= xx(i) <= xx(i3)=x3
1212 c-----------------------------------------------------------------------
1215 1 xx(i)=x0+(i-1.)/(i1-1.)*(x1-x0)
1217 2 xx(i)=x1+(i-i1*1.)/(i2-i1*1.)*(x2-x1)
1219 3 xx(i)=x2+(i-i2*1.)/(i3-i2*1.)*(x3-x2)
1223 cc---------------------------------------------------------------------
1224 c subroutine utaxis(i,j,a1,a2,a3)
1225 cc-----------------------------------------------------------------------
1226 cc calculates the axis defined by the ptls i,j in the i,j cm system
1227 cc---------------------------------------------------------------------
1228 c include 'epos.inc'
1229 c double precision pi1,pi2,pi3,pi4,pj1,pj2,pj3,pj4,p1,p2,p3,p4,p5
1234 c pi1=dble(pptl(1,i))
1235 c pi2=dble(pptl(2,i))
1236 c pi3=dble(pptl(3,i))
1237 c pi4=dble(pptl(4,i))
1238 c pj1=dble(pptl(1,j))
1239 c pj2=dble(pptl(2,j))
1240 c pj3=dble(pptl(3,j))
1241 c pj4=dble(pptl(4,j))
1246 c p5=dsqrt(p4**2-p3**2-p2**2-p1**2)
1247 c call utlob2(1,p1,p2,p3,p4,p5,pi1,pi2,pi3,pi4,50)
1248 c call utlob2(1,p1,p2,p3,p4,p5,pj1,pj2,pj3,pj4,51)
1249 c err=(pi1+pj1)**2+(pi2+pj2)**2+(pi3+pj3)**2
1250 c if(err.gt.1d-3)then
1251 c call utmsg('utaxis')
1252 c write(ifch,*)'***** err=',err
1253 c write(ifch,*)'pi:',pi1,pi2,pi3,pi4
1254 c write(ifch,*)'pj:',pj1,pj2,pj3,pj4
1257 c a=dsqrt( (pj1-pi1)**2 + (pj2-pi2)**2 + (pj3-pi3)**2 )
1258 c if(a.eq.0.d0)return
1259 c a1=sngl((pi1-pj1)/a)
1260 c a2=sngl((pi2-pj2)/a)
1261 c a3=sngl((pi3-pj3)/a)
1265 cc-----------------------------------------------------------------------
1266 c subroutine utchm(arp,arm,ii)
1267 cc-----------------------------------------------------------------------
1268 cc checks whether arp**2=0 and arm**2=0.
1269 cc-----------------------------------------------------------------------
1270 c include 'epos.inc'
1271 c double precision arp(4),arm(4),difp,difm
1272 c difp=arp(4)**2-arp(1)**2-arp(2)**2-arp(3)**2
1273 c difm=arm(4)**2-arm(1)**2-arm(2)**2-arm(3)**2
1274 c if(dabs(difp).gt.1e-3*arp(4)**2
1275 c *.or.dabs(difm).gt.1e-3*arm(4)**2)then
1276 c call utmsg('utchm ')
1277 c write(ifch,*)'***** mass non zero - ',ii
1278 c write(ifch,*)'jet-mass**2`s: ',difp,difm
1279 c write(ifch,*)'energy**2`s: ',arp(4)**2,arm(4)**2
1280 c write(ifch,*)(sngl(arp(i)),i=1,4)
1281 c write(ifch,*)(sngl(arm(i)),i=1,4)
1287 c-----------------------------------------------------------------------
1288 subroutine utclea(nptlii,nptl0)
1289 c-----------------------------------------------------------------------
1290 c starting from nptlii
1291 c overwrites istptl=99 particles in /cptl/, reduces so nptl
1292 c and update minfra and maxfra
1293 c-----------------------------------------------------------------------
1295 integer newptl(mxptl)!,oldptl(mxptl),ii(mxptl)
1298 if(ishsub/100.eq.18)ish=mod(ishsub,100)
1300 call utpri('utclea',ish,ishini,2)
1302 nptli=max(maproj+matarg+1,nptlii)
1307 if(ish.ge.2)write(ifch,*)'entering subr utclea:',nptl
1310 write(ifch,*)('-',l=1,68)
1311 write(ifch,*)'sr utclea. initial.'
1312 write(ifch,*)('-',l=1,68)
1314 write(ifch,116)iorptl(n),jorptl(n),n,ifrptl(1,n),ifrptl(2,n)
1315 *,idptl(n),sqrt(pptl(1,n)**2+pptl(2,n)**2),pptl(3,n),pptl(5,n)
1316 *,istptl(n),ityptl(n)
1318 116 format(1x,i6,i6,4x,i6,4x,i6,i6,i12,3x,3(e8.2,1x),i3,i3)
1323 c if(ishsub/100.eq.18)ish=mod(ishsub,100)
1327 if(i.gt.nptl)goto 1000
1328 if(istptl(i).eq.99)goto 2
1339 if(istptl(j).eq.99)goto 4
1342 c write(ifch,*)'move',j,' to ',i
1343 c write(ifch,*)idptl(i),ityptl(i),idptl(j),ityptl(j),minfra,maxfra
1345 if(j.ge.minfra0.and.j.le.maxfra0)then
1346 minfra1=min(minfra1,i)
1347 maxfra1=max(maxfra1,i)
1362 if(nptl0.gt.0)goto 20
1368 c if(io.le.0)ii(k)=io
1369 c if(io.gt.0)ii(k)=newptl(io)
1376 c if(jo.le.0)ii(k)=jo
1377 c if(jo.gt.0)ii(k)=newptl(jo)
1384 c if(if1.le.0)ii(k)=if1
1385 c if(if1.gt.0)ii(k)=newptl(if1)
1388 c16 ifrptl(1,k)=ii(k)
1392 c if(if2.le.0)ii(k)=if2
1393 c if(if2.gt.0)ii(k)=newptl(if2)
1396 c18 ifrptl(2,k)=ii(k)
1399 c if(ifrptl(1,k).eq.0.and.ifrptl(2,k).gt.0)ifrptl(1,k)=ifrptl(2,k)
1400 c if(ifrptl(2,k).eq.0.and.ifrptl(1,k).gt.0)ifrptl(2,k)=ifrptl(1,k)
1405 if(minfra1.lt.minfra0)minfra=minfra1
1406 if(maxfra1.ge.minfra1)maxfra=maxfra1
1409 write(ifch,*)'before exiting subr utclea:'
1411 write(ifch,116)iorptl(n),jorptl(n),n,ifrptl(1,n),ifrptl(2,n)
1412 *,idptl(n),sqrt(pptl(1,n)**2+pptl(2,n)**2),pptl(3,n),pptl(5,n)
1413 *,istptl(n),ityptl(n)
1417 if(ish.ge.2)write(ifch,*)'exiting subr utclea:',nptl
1420 call utprix('utclea',ish,ishini,2)
1425 c---------------------------------------------------------------------
1426 subroutine utfit(x,y,ndata,sig,mwt,a,b,siga,sigb,chi2,q)
1427 c---------------------------------------------------------------------
1428 c linear fit to data
1430 c ndata: nr of data points
1431 c x(),y(),sig(): data
1432 c mwt: unweighted (0) or weighted (else) data points
1434 c a,b: parameters of linear fit a+b*x
1435 c---------------------------------------------------------------------
1437 REAL a,b,chi2,q,siga,sigb,sig(ndata),x(ndata),y(ndata)
1440 REAL sigdat,ss,st2,sx,sxoss,sy,t,wt,utgmq
1463 t=(x(i)-sxoss)/sig(i)
1476 siga=sqrt((1.+sx*sx/(ss*st2))/ss)
1481 chi2=chi2+(y(i)-a-b*x(i))**2
1484 sigdat=sqrt(chi2/(ndata-2))
1489 chi2=chi2+((y(i)-a-b*x(i))/sig(i))**2
1491 q=utgmq(0.5*(ndata-2),0.5*chi2)
1496 c-----------------------------------------------------------------------
1498 c-----------------------------------------------------------------------
1499 c gamma fctn tabulated
1501 c-----------------------------------------------------------------------
1502 double precision utgamtab,utgam,al,dl
1503 common/gamtab/utgamtab(10000)
1505 if(x.gt.0.01.and.x.lt.99.99)then
1510 utgam1=real(utgamtab(k2)*dl+utgamtab(k1)*(1.d0-dl))
1514 utgam1=real(utgam(dble(x)))
1519 c-----------------------------------------------------------------------
1520 double precision function utgam2(x)
1521 c-----------------------------------------------------------------------
1522 c gamma fctn tabulated
1524 c-----------------------------------------------------------------------
1525 double precision utgamtab,x,al,dl,utgam
1526 common/gamtab/utgamtab(10000)
1528 if(x.gt.0.01d0.and.x.le.99.99d0)then
1533 utgam2=utgamtab(k2)*dl+utgamtab(k1)*(1.d0-dl)
1534 elseif(x.eq.0.d0)then
1542 c-----------------------------------------------------------------------
1543 double precision function utgam(x)
1544 c-----------------------------------------------------------------------
1547 c-----------------------------------------------------------------------
1549 double precision c(13),x,z,f
1551 1/ 0.00053 96989 58808, 0.00261 93072 82746, 0.02044 96308 23590,
1552 2 0.07309 48364 14370, 0.27964 36915 78538, 0.55338 76923 85769,
1553 3 0.99999 99999 99998,-0.00083 27247 08684, 0.00469 86580 79622,
1554 4 0.02252 38347 47260,-0.17044 79328 74746,-0.05681 03350 86194,
1555 5 1.13060 33572 86556/
1558 if(x .gt. 170.d0) goto6
1559 if(x .gt. 0.0d0) goto1
1560 if(x .eq. int(x)) goto5
1563 if(z .le. 1.0d0) goto4
1566 if(z .lt. 2.0d0) goto3
1572 1 f*((((((c(1)*z+c(2))*z+c(3))*z+c(4))*z+c(5))*z+c(6))*z+c(7))/
1573 2 ((((((c(8)*z+c(9))*z+c(10))*z+c(11))*z+c(12))*z+c(13))*z+1.0d0)
1574 if(x .gt. 0.0d0) return
1575 utgam=3.141592653589793d0/(sin(3.141592653589793d0*x)*utgam)
1577 5 write(ifch,10)sngl(x)
1578 10 format(1x,'argument of gamma fctn = ',e20.5)
1579 call utstop('utgam : negative integer argument&')
1580 6 write(ifch,11)sngl(x)
1581 11 format(1x,'argument of gamma fctn = ',e20.5)
1582 call utstop('utgam : argument too large&')
1585 c---------------------------------------------------------------------
1586 subroutine utgcf(gammcf,a,x,gln)
1587 c---------------------------------------------------------------------
1589 REAL a,gammcf,gln,x,EPS,FPMIN
1590 PARAMETER (ITMAX=100,EPS=3.e-7,FPMIN=1.e-30)
1593 REAL an,b,c,d,del,h,utgmln
1603 if(abs(d).lt.FPMIN)d=FPMIN
1605 if(abs(c).lt.FPMIN)c=FPMIN
1609 if(abs(del-1.).lt.EPS)goto 1
1611 pause 'a too large, ITMAX too small in utgcf'
1612 1 gammcf=exp(-x+a*log(x)-gln)*h
1616 c---------------------------------------------------------------------
1618 c---------------------------------------------------------------------
1621 DOUBLE PRECISION ser,stp,tmp,x,y,cof(6)
1623 DATA cof,stp/76.18009172947146d0,-86.50532032941677d0,
1624 *24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,
1625 *-.5395239384953d-5,2.5066282746310005d0/
1629 tmp=(x+0.5d0)*log(tmp)-tmp
1630 ser=1.000000000190015d0
1635 utgmln=tmp+log(stp*ser/x)
1639 c---------------------------------------------------------------------
1641 c---------------------------------------------------------------------
1643 CU USES utgcf,utgser
1644 REAL gammcf,gamser,gln
1645 if(x.lt.0..or.a.le.0.)pause 'bad arguments in utgmq'
1647 call utgser(gamser,a,x,gln)
1650 call utgcf(gammcf,a,x,gln)
1656 c---------------------------------------------------------------------
1657 subroutine utgser(gamser,a,x,gln)
1658 c---------------------------------------------------------------------
1660 REAL a,gamser,gln,x,EPS
1661 PARAMETER (ITMAX=100,EPS=3.e-7)
1664 REAL ap,del,sum,utgmln
1667 if(x.lt.0.)pause 'x < 0 in utgser'
1678 if(abs(del).lt.abs(sum)*EPS)goto 1
1680 pause 'a too large, ITMAX too small in utgser'
1681 1 gamser=sum*exp(-x+a*log(x)-gln)
1685 c-------------------------------------------------------------------------
1686 subroutine uticpl(ic,ifla,iqaq,iret)
1687 c-------------------------------------------------------------------------
1688 c adds a quark (iqaq=1) or antiquark (iqaq=2) of flavour ifla
1690 c-------------------------------------------------------------------------
1692 integer jc(nflav,2),ic(2)
1696 if(ish.ge.8)write(ifch,'(2i8,12i3)')ic,jc
1698 if(jc(ifla,jqaq).gt.0)then
1699 jc(ifla,jqaq)=jc(ifla,jqaq)-1
1701 jc(ifla,iqaq)=jc(ifla,iqaq)+1
1704 call idenco(jc,ic,ireten)
1705 if(ish.ge.8)write(ifch,'(2i8,12i3)')ic,jc
1706 if(ireten.eq.1)iret=1
1707 if(ic(1).eq.0.and.ic(2).eq.0.and.ireten.eq.0)then
1714 cc-----------------------------------------------------------------------
1715 c subroutine utindx(n,xar,x,i)
1716 cc-----------------------------------------------------------------------
1717 cc input: dimension n
1718 cc array xar(n) with xar(i) > xar(i-1)
1719 cc some number x between xar(1) and xar(n)
1720 cc output: the index i such that x is between xar(i) and xar(i+1)
1721 cc-----------------------------------------------------------------------
1722 c include 'epos.inc'
1724 c if(x.lt.xar(1))then
1726 c call utmsg('utindx')
1727 c write(ifch,*)'***** x=',x,' < xar(1)=',xar(1)
1732 c elseif(x.gt.xar(n))then
1734 c call utmsg('utindx')
1735 c write(ifch,*)'***** x=',x,' > xar(n)=',xar(n)
1744 c if((xar(lu).le.x).and.(x.le.xar(lz)))then
1746 c elseif((xar(lz).lt.x).and.(x.le.xar(lo)))then
1749 c call utstop('utindx: no interval found&')
1751 c if((lo-lu).ge.2) goto1
1752 c if(lo.le.lu)call utstop('utinvt: lo.le.lu&')
1757 c-----------------------------------------------------------------------
1758 function utinvt(n,x,q,y)
1759 c-----------------------------------------------------------------------
1760 c returns x with y=q(x)
1761 c-----------------------------------------------------------------------
1764 if(q(n).eq.0.)call utstop('utinvt: q(n)=0&')
1767 call utmsg('utinvt')
1768 write(ifch,*)'***** y=',y,' < 0'
1772 elseif(y.gt.q(n))then
1774 call utmsg('utinvt')
1775 write(ifch,*)'***** y=',y,' > ',q(n)
1783 if((q(lu).le.y).and.(y.le.q(lz)))then
1785 elseif((q(lz).lt.y).and.(y.le.q(lo)))then
1788 write(ifch,*)'q(1),y,q(n):',q(1),y,q(n)
1789 write(ifch,*)'lu,lz,lo:',lu,lz,lo
1790 write(ifch,*)'q(lu),q(lz),q(lo):',q(lu),q(lz),q(lo)
1791 call utstop('utinvt: no interval found&')
1793 if((lo-lu).ge.2) goto1
1794 if(lo.le.lu)call utstop('utinvt: lo.le.lu&')
1795 utinvt=x(lu)+(y-q(lu))*(x(lo)-x(lu))/(q(lo)-q(lu))
1799 c-----------------------------------------------------------------------
1800 subroutine utlob2(isig,p1,p2,p3,p4,p5,x1,x2,x3,x4,idi)
1801 c-----------------------------------------------------------------------
1802 c performs a lorentz boost, double prec.
1803 c isig=+1 is to boost the four vector x1,x2,x3,x4 such as to obtain it
1804 c in the frame specified by the 5-vector p1...p5 (5-vector=4-vector+mass).
1805 c isig=-1: the other way round, that means,
1806 c if the 4-vector x1...x4 is given in some frame characterized by
1807 c p1...p5 with respect to to some lab-frame, utlob2 returns the 4-vector
1808 c x1...x4 in the lab frame.
1809 c idi is a call identifyer (integer) to identify the call in case of problem
1810 c-----------------------------------------------------------------------
1812 double precision beta(4),z(4),p1,p2,p3,p4,p5,pp,bp,x1,x2,x3,x4
1813 *,xx0,x10,x20,x30,x40,x4x,x0123
1816 write(ifch,101)x1,x2,x3,x4,(x4-x3)*(x4+x3)-x2*x2-x1*x1
1817 write(ifch,301)p1,p2,p3,p4,p5,(p4-p3)*(p4+p3)-p2*p2-p1*p1
1818 101 format(' utlob2: x = ',5e13.5)
1819 301 format(' p = ',6e13.5)
1821 pp=(p4-p3)*(p4+p3)-p2*p2-p1*p1
1822 if(dabs(pp-p5*p5).gt.1e-3*p4*p4.and.dabs(pp-p5*p5).gt.1e-3)then
1823 call utmsg('utlob2')
1824 write(ifch,*)'***** p**2 .ne. p5**2'
1825 write(ifch,*)'call identifyer:',idi
1826 write(ifch,*)'p**2,p5**2: ',pp,p5*p5
1827 write(ifch,*)'p: ',p1,p2,p3,p4,p5
1835 xx0=(x4-x3)*(x4+x3)-x2*x2-x1*x1
1837 call utmsg('utlob2')
1838 write(ifch,*)'***** p5 negative.'
1839 write(ifch,*)'call identifyer:',idi
1840 write(ifch,*)'p(5): ',p1,p2,p3,p4,p5
1841 write(ifmt,*)'call identifyer:',idi
1842 write(ifmt,*)'p(5): ',p1,p2,p3,p4,p5
1844 call utstop('utlob2: p5 negative.&')
1856 220 bp=bp+z(k)*isig*beta(k)
1858 230 z(k)=z(k)+isig*beta(k)*z(4)
1859 *+isig*beta(k)*bp/(beta(4)+1.)
1860 z(4)=beta(4)*z(4)+bp
1866 *write(ifch,101)x1,x2,x3,x4,(x4-x3)*(x4+x3)-x2*x2-x1*x1
1868 x0123=xx0+x1*x1+x2*x2+x3*x3
1870 x4=sign( dsqrt(x0123) , x4x )
1875 write(ifch,101)x1,x2,x3,x4,(x4-x3)*(x4+x3)-x2*x2-x1*x1
1878 if(ish.ge.9)write(ifch,*)'check x**2_ini -- x**2_fin'
1879 if(dabs(x4-x4x).gt.1d-2*dabs(x4).and.dabs(x4-x4x).gt.1d-2)then
1880 call utmsg('utlob2')
1881 write(ifch,*)'***** x**2_ini .ne. x**2_fin.'
1882 write(ifch,*)'call identifyer:',idi
1883 write(ifch,*)'x1 x2 x3 x4 x**2 (initial/final/corrected):'
1885 write(ifch,102)x10,x20,x30,x40,(x40-x30)*(x40+x30)-x20*x20-x10*x10
1886 write(ifch,102)x1,x2,x3,x4x,(x4x-x3)*(x4x+x3)-x2*x2-x1*x1
1887 write(ifch,102)x1,x2,x3,x4,(x4-x3)*(x4+x3)-x2*x2-x1*x1
1891 if(ish.ge.9)write(ifch,*)'return from utlob2'
1895 c-----------------------------------------------------------------------
1896 subroutine utlob3(isig,p1,p2,p3,p4,p5,x1,x2,x3,x4)
1897 c-----------------------------------------------------------------------
1898 c performs a lorentz boost, double prec.
1899 c but arguments are single precision
1900 c-----------------------------------------------------------------------
1901 double precision xx1,xx2,xx3,xx4
1907 *,dble(p1),dble(p2),dble(p3),dble(p4),dble(p5)
1908 *,xx1,xx2,xx3,xx4,52)
1916 c-----------------------------------------------------------------------
1917 subroutine utlob5(yboost,x1,x2,x3,x4,x5)
1918 c-----------------------------------------------------------------------
1919 amt=sqrt(x5**2+x1**2+x2**2)
1920 y=sign(1.,x3)*alog((x4+abs(x3))/amt)
1927 c-----------------------------------------------------------------------
1928 subroutine utlob4(isig,pp1,pp2,pp3,pp4,pp5,x1,x2,x3,x4)
1929 c-----------------------------------------------------------------------
1930 c performs a lorentz boost, double prec.
1931 c but arguments are partly single precision
1932 c-----------------------------------------------------------------------
1933 double precision xx1,xx2,xx3,xx4,pp1,pp2,pp3,pp4,pp5
1938 call utlob2(isig,pp1,pp2,pp3,pp4,pp5,xx1,xx2,xx3,xx4,53)
1947 c-----------------------------------------------------------------------
1948 subroutine utlobo(isig,p1,p2,p3,p4,p5,x1,x2,x3,x4)
1949 c-----------------------------------------------------------------------
1950 c performs a lorentz boost
1951 c-----------------------------------------------------------------------
1955 call utmsg('utlobo')
1956 write(ifch,*)'***** mass <= 0.'
1957 write(ifch,*)'p(5): ',p1,p2,p3,p4,p5
1959 call utstop('utlobo: mass <= 0.&')
1971 220 bp=bp+z(k)*isig*beta(k)
1973 230 z(k)=z(k)+isig*beta(k)*z(4)
1974 *+isig*beta(k)*bp/(beta(4)+1.)
1975 z(4)=beta(4)*z(4)+bp
1983 c-----------------------------------------------------------------------
1984 subroutine utloc(ar,n,a,l)
1985 c-----------------------------------------------------------------------
1989 if(a.lt.ar(i))return
1995 cc-----------------------------------------------------------------------
1996 c subroutine utlow(cone)
1997 cc-----------------------------------------------------------------------
1999 c if(cone.eq.'A')cone='a'
2000 c if(cone.eq.'B')cone='b'
2001 c if(cone.eq.'C')cone='c'
2002 c if(cone.eq.'D')cone='d'
2003 c if(cone.eq.'E')cone='e'
2004 c if(cone.eq.'F')cone='f'
2005 c if(cone.eq.'G')cone='g'
2006 c if(cone.eq.'H')cone='h'
2007 c if(cone.eq.'I')cone='i'
2008 c if(cone.eq.'J')cone='j'
2009 c if(cone.eq.'K')cone='k'
2010 c if(cone.eq.'L')cone='l'
2011 c if(cone.eq.'M')cone='m'
2012 c if(cone.eq.'N')cone='n'
2013 c if(cone.eq.'O')cone='o'
2014 c if(cone.eq.'P')cone='p'
2015 c if(cone.eq.'Q')cone='q'
2016 c if(cone.eq.'R')cone='r'
2017 c if(cone.eq.'S')cone='s'
2018 c if(cone.eq.'T')cone='t'
2019 c if(cone.eq.'U')cone='u'
2020 c if(cone.eq.'V')cone='v'
2021 c if(cone.eq.'W')cone='w'
2022 c if(cone.eq.'X')cone='x'
2023 c if(cone.eq.'Y')cone='y'
2024 c if(cone.eq.'Z')cone='z'
2028 cc-----------------------------------------------------------------------
2029 c subroutine utlow3(cthree)
2030 cc-----------------------------------------------------------------------
2031 c character cthree*3
2033 c1 call utlow(cthree(i:i))
2037 cc-----------------------------------------------------------------------
2038 c subroutine utlow6(csix)
2039 cc-----------------------------------------------------------------------
2042 c1 call utlow(csix(i:i))
2046 cc-----------------------------------------------------------------------
2047 c function utmom(k,n,x,q)
2048 cc-----------------------------------------------------------------------
2049 cc calculates kth moment for f(x) with q(i)=int[0,x(i)]f(z)dz
2050 cc-----------------------------------------------------------------------
2052 c if(n.lt.2)call utstop('utmom : dimension too small&')
2055 c1 utmom=utmom+((x(i)+x(i-1))/2)**k*(q(i)-q(i-1))
2060 c-----------------------------------------------------------------------
2061 function utpcm(a,b,c)
2062 c-----------------------------------------------------------------------
2063 c calculates cm momentum for a-->b+c
2064 c-----------------------------------------------------------------------
2065 val=(a*a-b*b-c*c)*(a*a-b*b-c*c)-(2.*b*c)*(2.*b*c)
2066 if(val.lt.0..and.val.gt.-1e-4)then
2070 utpcm=sqrt(val)/(2.*a)
2074 c-----------------------------------------------------------------------
2075 double precision function utpcmd(a,b,c)
2076 c-----------------------------------------------------------------------
2077 c calculates cm momentum for a-->b+c
2078 c-----------------------------------------------------------------------
2079 double precision a,b,c,val
2080 val=(a*a-b*b-c*c)*(a*a-b*b-c*c)-(2.*b*c)*(2.*b*c)
2081 if(val.lt.0d0.and.val.gt.-1d-4)then
2084 elseif(val.lt.0d0)then
2085 call utstop('utpcmd: problem with masses&')
2087 utpcmd=sqrt(val)/(2.d0*a)
2091 c-----------------------------------------------------------------------
2092 subroutine utpri(text,idmmy,ishini,ishx)
2093 c-----------------------------------------------------------------------
2096 c double precision seedx !!!
2098 if(ishevt.ne.0.and.nrevt+1.ne.ishevt)return
2101 if(subpri(nr)(1:6).eq.text)then
2107 write(ifch,'(1x,43a)')
2108 * ('-',i=1,10),' entry ',text,' ',('-',i=1,30)
2109 c call ranfgt(seedx) !!!
2110 c if(ish.ge.ishx)write(ifch,*)'seed:',seedx !!!
2115 c-----------------------------------------------------------------------
2116 subroutine utprix(text,idmmy,ishini,ishx)
2117 c-----------------------------------------------------------------------
2120 if(ishevt.ne.0.and.nrevt+1.ne.ishevt)return
2121 if(ish.ge.ishx)write(ifch,'(1x,44a)')
2122 *('-',i=1,30),' exit ',text,' ',('-',i=1,11)
2127 c-----------------------------------------------------------------------
2128 subroutine utprj(text,idmmy,ishini,ishx)
2129 c-----------------------------------------------------------------------
2132 c double precision seedx !!!
2133 idx=index(text,' ')-1
2135 if(ishevt.ne.0.and.nrevt+1.ne.ishevt)return
2138 if(subpri(nr)(1:idx).eq.text(1:idx))then
2144 write(ifch,'(1x,43a)')
2145 * ('-',i=1,10),' entry ',text(1:idx),' ',('-',i=1,30)
2146 c call ranfgt(seedx) !!!
2147 c if(ish.ge.ishx)write(ifch,*)'seed:',seedx !!!
2152 c-----------------------------------------------------------------------
2153 subroutine utprjx(text,idmmy,ishini,ishx)
2154 c-----------------------------------------------------------------------
2157 idx=index(text,' ')-1
2158 if(ishevt.ne.0.and.nrevt+1.ne.ishevt)return
2159 if(ish.ge.ishx)write(ifch,'(1x,44a)')
2160 *('-',i=1,30),' exit ',text(1:idx),' ',('-',i=1,11)
2165 c-----------------------------------------------------------------------
2166 function utquad(m,x,f,k)
2167 c-----------------------------------------------------------------------
2168 c performs an integration according to simpson
2169 c-----------------------------------------------------------------------
2173 1 utquad=utquad+(f(i)+f(i+1))/2*(x(i+1)-x(i))
2177 c-----------------------------------------------------------------------
2178 subroutine utquaf(fu,n,x,q,wo,x0,x1,x2,x3)
2179 c-----------------------------------------------------------------------
2180 c returns q(i) = integral [x(1)->x(i)] fu(x) dx
2181 c-----------------------------------------------------------------------
2183 real x(n),q(n),wo(n)
2187 if(x1.lt.x0.or.x2.lt.x1.or.x3.lt.x2)then
2189 call utmsg('utquaf')
2190 write(ifch,*)' xi=',x0,x1,x2,x3
2194 call utar(n/3,n*2/3,n,x0,x1,x2,x3,x)
2198 z=x(i-1)+(k-1.)/(m-1.)*(x(i)-x(i-1))
2201 q(i)=q(i-1)+utquad(m,xa,fa,m)
2206 c-----------------------------------------------------------------------
2207 subroutine utrepl(i,j)
2208 c-----------------------------------------------------------------------
2209 c i is replaced by j in /cptl/
2210 c-----------------------------------------------------------------------
2213 1 pptl(k,i) =pptl(k,j)
2214 iorptl(i) = 0 !iorptl(j)
2216 istptl(i) =istptl(j)
2218 2 tivptl(k,i)=tivptl(k,j)
2220 3 ifrptl(k,i)= 0 !ifrptl(k,j)
2221 jorptl(i) = 0 !jorptl(j)
2223 4 xorptl(k,i)=xorptl(k,j)
2225 5 ibptl(k,i) =ibptl(k,j)
2226 ityptl(i) =ityptl(j)
2227 iaaptl(i) =iaaptl(j)
2228 radptl(i) =radptl(j)
2229 desptl(i) =desptl(j)
2230 dezptl(i) =dezptl(j)
2231 qsqptl(i) =qsqptl(j)
2232 zpaptl(i) =zpaptl(j)
2233 itsptl(i) =itsptl(j)
2234 rinptl(i) =rinptl(j)
2238 cc-----------------------------------------------------------------------
2239 c subroutine utresm(icp1,icp2,icm1,icm2,amp,idpr,iadj,ireten)
2240 cc-----------------------------------------------------------------------
2241 c parameter (nflav=6)
2242 c integer icm(2),icp(2),jcm(nflav,2),jcp(nflav,2)
2247 c CALL IDDECO(ICM,JCM)
2248 c CALL IDDECO(ICP,JCP)
2251 c37 jcP(nf,k)=jcp(nf,k)+jcm(nf,k)
2252 c CALL IDENCO(JCP,ICP,IRETEN)
2253 c IDP=IDTRA(ICP,0,0,3)
2254 c call idres(idp,amp,idpr,iadj,0)
2258 cc-----------------------------------------------------------------------
2259 c subroutine utroa1(phi,a1,a2,a3,x1,x2,x3)
2260 cc-----------------------------------------------------------------------
2261 cc rotates x by angle phi around axis a.
2262 cc normalization of a is irrelevant.
2263 cc-----------------------------------------------------------------------
2264 c double precision aaa,aa(3),xxx,xx(3),e1(3),e2(3),e3(3),xp,xt,dphi
2278 c if(xxx.eq.0d0)return
2279 c if(aaa.eq.0d0)call utstop('utroa1: zero rotation axis&')
2292 c if(xxx**2-xp**2.le.0.)return
2293 c xt=dsqrt(xxx**2-xp**2)
2294 cc e1 = vector x_transverse / absolute value x_transverse
2296 c e1(i)=(xx(i)-e3(i)*xp)/xt
2298 cc e2 orthogonal e3,e1
2299 c call utvec2(e3,e1,e2)
2302 c xx(i)=xp*e3(i)+xt*dcos(dphi)*e1(i)+xt*dsin(dphi)*e2(i)
2304 cc back to single precision
2311 cc-----------------------------------------------------------------------
2312 c subroutine utroa2(phi,a1,a2,a3,x1,x2,x3)
2313 cc-----------------------------------------------------------------------
2314 cc rotates x by angle phi around axis a.
2315 cc normalization of a is irrelevant.
2316 cc double precision phi,a1,a2,a3,x1,x2,x3
2317 cc-----------------------------------------------------------------------
2318 c double precision phi,a1,a2,a3,x1,x2,x3
2319 c double precision aaa,aa(3),xxx,xx(3),e1(3),e2(3),e3(3),xp,xt,dphi
2333 c if(xxx.eq.0d0)return
2334 c if(aaa.eq.0d0)call utstop('utroa1: zero rotation axis&')
2335 c aaa=1.0/dsqrt(aaa)
2346 c if(xxx-xp**2.le.0.)return
2347 c xt=dsqrt(xxx-xp**2)
2348 cc e1 = vector x_transverse / absolute value x_transverse
2350 c e1(i)=(xx(i)-e3(i)*xp)/xt
2352 cc e2 orthogonal e3,e1
2353 c call utvec2(e3,e1,e2)
2356 c xx(i)=xp*e3(i)+xt*dcos(dphi)*e1(i)+xt*dsin(dphi)*e2(i)
2358 cc back to single precision
2365 cc--------------------------------------------------------------------
2366 c function utroot(funcd,x1,x2,xacc)
2367 cc--------------------------------------------------------------------
2368 cc combination of newton-raphson and bisection method for root finding
2370 cc funcd: subr returning fctn value and first derivative
2371 cc x1,x2: x-interval
2375 cc--------------------------------------------------------------------
2376 c include 'epos.inc'
2378 c REAL utroot,x1,x2,xacc
2380 c PARAMETER (MAXIT=100)
2382 c REAL df,dx,dxold,f,fh,fl,temp,xh,xl
2383 c call funcd(x1,fl,df)
2384 c call funcd(x2,fh,df)
2385 c if((fl.gt.0..and.fh.gt.0.).or.(fl.lt.0..and.fh.lt.0.))
2386 c *call utstop('utroot: root must be bracketed&')
2390 c else if(fh.eq.0.)then
2393 c else if(fl.lt.0.)then
2403 c call funcd(utroot,f,df)
2405 c if(((utroot-xh)*df-f)*((utroot-xl)*df-f).ge.0..or. abs(2.*
2406 c *f).gt.abs(dxold*df) ) then
2410 c if(xl.eq.utroot)return
2416 c if(temp.eq.utroot)return
2418 c if(abs(dx).lt.xacc) return
2419 c call funcd(utroot,f,df)
2426 c call utmsg('utroot')
2427 c write(ifch,*)'***** exceeding maximum iterations'
2428 c write(ifch,*)'dx:',dx
2433 c-----------------------------------------------------------------------
2434 subroutine utrot2(isig,ax,ay,az,x,y,z)
2435 c-----------------------------------------------------------------------
2436 c performs a rotation, double prec.
2437 c-----------------------------------------------------------------------
2439 double precision ax,ay,az,x,y,z,rx,ry,rz
2440 *,alp,bet,cosa,sina,cosb,sinb,xs,ys,zs
2441 if(ax**2.eq.0.and.ay**2.eq.0.and.az**2.eq.0.)then
2442 write(ifch,*)'ax**2,ay**2,az**2:',ax**2,ay**2,az**2
2443 write(ifch,*)'ax,ay,az:',ax,ay,az
2444 call utstop('utrot2: zero vector.&')
2455 if(rz**2+ry**2.ne.0.)then
2456 alp=dabs(dacos(rz/dsqrt(rz**2+ry**2)))*sign(1.,sngl(ry))
2458 *dabs(dacos(dsqrt(rz**2+ry**2)/dsqrt(rz**2+ry**2+rx**2)))*
2469 xs=x*cosb-y*sina*sinb-z*cosa*sinb
2471 zs=x*sinb+y*sina*cosb+z*cosa*cosb
2472 elseif(isig.lt.0)then
2474 ys=-x*sinb*sina+y*cosa+z*cosb*sina
2475 zs=-x*sinb*cosa-y*sina+z*cosb*cosa
2484 c-----------------------------------------------------------------------
2485 subroutine utrot4(isig,ax,ay,az,x,y,z)
2486 c-----------------------------------------------------------------------
2487 c performs a rotation, double prec.
2488 c arguments partly single
2489 c-----------------------------------------------------------------------
2490 double precision ax,ay,az,xx,yy,zz
2494 call utrot2(isig,ax,ay,az,xx,yy,zz)
2501 c-----------------------------------------------------------------------
2502 subroutine utrota(isig,ax,ay,az,x,y,z)
2503 c-----------------------------------------------------------------------
2504 c performs a rotation
2505 c-----------------------------------------------------------------------
2515 if(rz.eq.0..and.ry.eq.0.)then
2519 alp=abs(utacos(rz/sqrt(rz**2+ry**2)))*sign(1.,ry)
2522 *abs(utacos(sqrt(rz**2+ry**2)/sqrt(rz**2+ry**2+rx**2)))*sign(1.,rx)
2528 xs=x*cosb-y*sina*sinb-z*cosa*sinb
2530 zs=x*sinb+y*sina*cosb+z*cosa*cosb
2531 elseif(isig.lt.0)then
2533 ys=-x*sinb*sina+y*cosa+z*cosb*sina
2534 zs=-x*sinb*cosa-y*sina+z*cosb*cosa
2542 c-----------------------------------------------------------------------
2543 subroutine utstop(text)
2544 c-----------------------------------------------------------------------
2545 c returns error message and stops execution.
2546 c text is an optonal text to appear in the error message.
2547 c text is a character string of length 40;
2548 c for shorter text, it has to be terminated by &;
2549 c example: call utstop('error in subr xyz&')
2550 c-----------------------------------------------------------------------
2553 character text*40 ,txt*6
2556 if(text(i:i).eq.'&')imax=i
2561 write(ifi,101)('*',k=1,72),text(1:imax-1)
2562 *,nrevt+1,seedc,('*',k=1,72)
2565 */1x,'***** stop in ',a
2566 */1x,'***** current event number: ',i12
2567 */1x,'***** initial seed for current event:',d25.15
2573 write(ifch,'(1x,74a1)')('*',j=1,72)
2574 write(ifch,100)txt,nrevt+1,seedc
2575 100 format(1x,'***** msg from ',a6,'. es:',i7,2x,d23.17)
2579 write(ifch,'(1x,74a1)')('*',j=1,72)
2582 c-----------------------------------------------------------------
2583 subroutine uttrap(func,a,b,s)
2584 c-----------------------------------------------------------------
2585 c trapezoidal method for integration.
2586 c input: fctn func and limits a,b
2587 c output: value s of the integral
2588 c-----------------------------------------------------------------
2592 REAL a,b,func,s,epsr
2600 if(ish.ge.9)write(ifch,*)'sr uttrap: j:',j
2601 call uttras(func,a,b,s,j)
2603 if (ds.lt.epsr*abs(olds)) return
2609 call utmsg('uttrap')
2611 *'***** requested accuracy could not be achieved'
2612 write(ifch,*)'achieved accuracy: ',ds/abs(olds)
2613 write(ifch,*)'requested accuracy:',epsr
2619 c-----------------------------------------------------------------
2620 subroutine uttraq(func,a,b,s)
2621 c-----------------------------------------------------------------
2622 c trapezoidal method for integration.
2623 c input: function func and limits a,b
2624 c output: value s of the integral
2625 c-----------------------------------------------------------------
2629 PARAMETER (eps=1.e-6)
2635 10 call uttras(func,a,b,s,j)
2637 if (ds.le.eps*abs(olds)) return
2644 c-----------------------------------------------------------------
2645 subroutine uttras(func,a,b,s,n)
2646 c-----------------------------------------------------------------
2647 c performs one iteration of the trapezoidal method for integration
2648 c-----------------------------------------------------------------
2655 s=0.5*(b-a)*(func(a)+func(b))
2666 s=0.5*(s+(b-a)*sum/tnm)
2671 cc-----------------------------------------------------------------------
2672 c subroutine utvec1(a,b,c)
2673 cc-----------------------------------------------------------------------
2674 cc returns vector product c = a x b .
2675 cc-----------------------------------------------------------------------
2676 c real a(3),b(3),c(3)
2677 c c(1)=a(2)*b(3)-a(3)*b(2)
2678 c c(2)=a(3)*b(1)-a(1)*b(3)
2679 c c(3)=a(1)*b(2)-a(2)*b(1)
2683 cc-----------------------------------------------------------------------
2684 c subroutine utvec2(a,b,c)
2685 cc-----------------------------------------------------------------------
2686 cc returns vector product c = a x b .
2687 cc a,b,c double precision.
2688 cc-----------------------------------------------------------------------
2689 c double precision a(3),b(3),c(3)
2690 c c(1)=a(2)*b(3)-a(3)*b(2)
2691 c c(2)=a(3)*b(1)-a(1)*b(3)
2692 c c(3)=a(1)*b(2)-a(2)*b(1)
2696 c-------------------------------------------------------------------
2697 subroutine utword(line,i,j,iqu)
2698 c-------------------------------------------------------------------
2699 c finds the first word of the character string line(j+1:160).
2700 c the word is line(i:j) (with new i and j).
2701 c if j<0 or if no word found --> new line read.
2702 c a text between quotes "..." is considered one word;
2703 c stronger: a text between double quotes ""..."" is consid one word
2704 c stronger: a text between "{ and }" is considered one word
2705 c-------------------------------------------------------------------
2707 c line: character string (*160)
2708 c i: integer between 1 and 160
2709 c iqu: for iqu=1 a "?" is written to output before reading a line,
2710 c otherwise (iqu/=1) nothing is typed
2712 c i,j: left and right end of word (word=line(i:j))
2713 c-------------------------------------------------------------------
2716 character*1 empty(mempty),mk
2720 parameter(mxdefine=40)
2721 character w1define*100,w2define*100
2722 common/cdefine/ndefine,l1define(mxdefine),l2define(mxdefine)
2723 & ,w1define(mxdefine),w2define(mxdefine)
2731 if(iqu.eq.1.and.iprmpt.gt.0)write(ifmt,'(a)')'?'
2732 if(nopen.eq.0)ifopx=ifop
2733 if(nopen.gt.0)ifopx=20+nopen
2734 if(nopen.lt.0)ifopx=ifcp
2735 read(ifopx,'(a160)',end=9999)line
2736 if(iecho.eq.1.or.(nopen.ge.0.and.kcpopen.eq.1))then
2739 if(line(k:k).ne.' ')kmax=k
2742 if(nopen.ge.0.and.kcpopen.eq.1)
2743 & write(ifcp,'(a)')line(1:kmax)
2745 & write(ifmt,'(a)')line(1:kmax)
2750 if(line(i:i).eq.'!')goto5
2752 if(line(i:i).eq.empty(ne))goto1
2757 if(line(i:i).eq.'~')mk='~'
2758 if(line(i:i+1).eq.'"{')mrk='}"'
2759 if(line(i:i+1).eq.'""')mrk='""'
2760 if(mrk.ne.' ')goto10
2761 if(line(i:i).eq.'"')mk='"'
2766 if(line(j:j).eq.'!')goto7
2768 if(line(j:j).eq.empty(ne))goto7
2773 if(i.ge.160-1)stop'utword: make line shorter!!! '
2776 if(line(j:j).eq.mk)stop'utword: empty string!!! '
2779 if(line(j:j).eq.mk)then
2787 if(i.ge.160-3)stop'utword: make line shorter!!!! '
2790 if(line(j:j+1).eq.mrk)stop'utword: empty string!!!! '
2793 if(line(j:j+1).eq.mrk)then
2801 !--------#define---------------
2802 if(ndefine.gt.0)then
2807 if(line(i0:i0-1+l1).eq.w1define(ndf)(1:l1))then
2809 line(i0:i0-1+l1)=w2define(ndf)(1:l2)
2810 elseif(l2.lt.l1)then
2811 line(i0:i0+l2-1)=w2define(ndf)(1:l2)
2815 elseif(l2.gt.l1)then
2817 if(line(k:k).ne.' ')
2818 & stop'utword: no space for `define` replacement. '
2820 line(i0:i0+l2-1)=w2define(ndf)(1:l2)
2827 if(line(k:k).ne.' ')j0=j
2836 if(nopen.eq.0.and.iprmpt.eq.-1)iprmpt=1
2840 c--------------------------------------------------------------------
2841 subroutine utworn(line,j,ne)
2842 c--------------------------------------------------------------------
2843 c returns number ne of nonempty characters of line(j+1:160)
2844 c--------------------------------------------------------------------
2848 if(line(l:l).ne.' ')ne=ne+1
2853 c-----------------------------------------------------------------------
2854 subroutine getairmol(iz,ia)
2855 c-----------------------------------------------------------------------
2859 do while(r.gt.0.) ! choose air-molecule
2863 iz = nint(airznxs(i))
2864 ia = nint(airanxs(i))
2867 c----------------------------------------------------------------------
2869 subroutine factoriel
2871 c----------------------------------------------------------------------
2872 c tabulation of fctrl(n)=n!, facto(n)=1/n! and utgamtab(x) for x=0 to 50
2873 c----------------------------------------------------------------------
2874 include 'epos.incems'
2875 double precision utgamtab,utgam,x
2876 common/gamtab/utgamtab(10000)
2881 do i=1,min(npommx,nfctrl)
2882 fctrl(i)=fctrl(i-1)*dble(i)
2883 facto(i)=1.d0/fctrl(i)
2888 utgamtab(k)=utgam(x)
2894 c-----------------------------------------------------------------------
2895 subroutine fremnu(amin,ca,cb,ca0,cb0,ic1,ic2,ic3,ic4)
2896 c-----------------------------------------------------------------------
2897 common/hadr2/iomodl,idproj,idtarg,wexcit
2899 common/hadr1/pnll,ptq,exmass,cutmss,wproj,wtarg
2901 integer ic(2),jc(nflav,2)
2911 amin=utamnu(keu,ked,kes,kec,keb,ket,5) !???4=2mults, 5=1mult
2912 if(ca-ca0.eq.0..and.cb-cb0.eq.0..and.rangen().gt.wexcit)then
2946 c-----------------------------------------------------------------------
2947 function fremnux(ic)
2948 c-----------------------------------------------------------------------
2950 common/hadr1/pnll,ptq,exmass,cutmss,wproj,wtarg
2952 integer ic(2),jc(nflav,2)
2962 fremnux=utamnu(keu,ked,kes,kec,keb,ket,4) !+exmass !???4=2mults, 5=1mult
2966 c-----------------------------------------------------------------------
2967 function fremnux2(ic)
2968 c-----------------------------------------------------------------------
2970 common/hadr1/pnll,ptq,exmass,cutmss,wproj,wtarg
2972 integer ic(2),jc(nflav,2)
2982 fremnux2=utamnu(keu,ked,kes,kec,keb,ket,5) !+exmass !???4=2mults, 5=1mult
2986 c-----------------------------------------------------------------------
2987 subroutine fremnx(ammax,amin,sm,ic3,ic4,iret)
2988 c-----------------------------------------------------------------------
2992 if(ic3.eq.0.and.ic4.eq.0)then
2993 if(ammax.lt.amin**2)then
2999 c ammax1=min(ammax,(engy/4.)**2)
3001 if(ammax1.lt.amin**2)then
3006 sm=amin**2*(ammax1/amin**2)**rangen()
3008 sm=amin**2*(1.+((ammax1/amin**2)**(1.+alpr)-1.)
3009 * *rangen())**(1./(1.+alpr))