c $Id: anndec.f,v 1.20 2003/05/02 13:14:47 weber Exp $
C####C##1#########2#########3#########4#########5#########6#########7##
subroutine anndec(io,mm1,ii1,iiz1,mm2,ii2,iiz2,sqrts,sig,gam)
c
c
cinput io : 0: annihilation; 1: decay
cinput mm1 : mass of scattering/decaying particle 1
cinput ii1 : ID of scattering/decaying particle 1
cinput iiz1 : $2\cdot I_3$ of scattering/decaying particle 1
cinput mm2 : mass of scattering particle 2
cinput ii2 : ID of scattering particle 2
cinput iiz2 : $2\cdot I_3$ of scattering particle 2
cinput sqrts : $sqrts{s}$ of collision; resonance mass for decay
coutput sig : resonance scattering cross section
coutput gam : width of the produced resonance
c
c {\tt anndec} handles meson-baryon and meson-meson annihilations
c as well as all meson and baryon resonance decays. In case of
c annihilations it returs the total resonance production cross section
c and the decay width of the resonance chosen as final state.
c The final state itself for both cases, annihilation and decay
c is returned via the {\tt newpart} common block. In the case
c of a decay the final state may consist of up to 4 particles.
c
c Technically, {\tt anndec} is only an interface to {\tt anndex},
c which actually handles the summation of the Breit-Wigner formulas
c in the annihilation case and the final state generation for
c annihilation and decay.
c
Ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C
implicit none
integer io,i1,i2,iz1,iz2,ii1,ii2,iiz1,iiz2,is
real*8 m1,m2,mm1,mm2,sig,gam,sqrts
include 'comres.f'
include 'options.f'
integer strit
C
C ************************************************************************
C Case 1 : Two ingoing Particles --> One outgoing Particle (Resonance,...)
C
C
i1=ii1
i2=ii2
iz1=iiz1
iz2=iiz2
m1=mm1
m2=mm2
if(io.eq.0)then !annihilation
C
sig=0d0
gam=0d0
C
C Check if (sqrt(s)-masses of ingoing particles) is significant different
C from zero
C
if(sqrts-mm1-mm2.le.1d-3)return
C
C Check if CTOption(15) is set different from zero-->then skip anndec.f
C
if(CTOption(15).ne.0)return
C
C
C Check if itype of particle one is smaller than the one of particle two
C if so --> interchange particle one and particle two
C in case of particle one = B and particle two = M --> then
C new particle one = M and new particle two = B
C
if(iabs(i1).lt.iabs(i2))call swpizm(i1,iz1,m1,i2,iz2,m2)
C
C Determination of the amount of the netstrangness
C
is=iabs(strit(i1)+strit(i2))
C
C maxbar (Maximum Baryon ityp)
C if second particle is antibaryon and first particle is strange
C switch antibaryon to baryon
C
if(iabs(i2).le.maxbar)then
if(i2.lt.0)then
if(strit(i1).ne.0)i1=-i1 ! get corresponding anti-branch
end if
end if
C
C
C Check if both particles are mesons
C
if(iabs(i1).ge.minmes.and.iabs(i2).ge.minmes)then
C
C
c... boson+boson sector
C
C
C Check if amount of netstrangeness is greater than 1
c currently no resonant processes for |s|>1 are implemented
if(is.gt.1)return
if(is.ne.0)then
call anndex(0,m1,i1,iz1,m2,i2,iz2,sqrts,
. sig,gam,maxbrm,minmes+1,maxmes,bmtype,branmes)
else
call anndex(0,m1,i1,iz1,m2,i2,iz2,sqrts,
. sig,gam,maxbrm,minmes+1,maxmes,bmtype,branmes)
endif
C
C Check if second particle is baryon
C (with zero amount of netstrangeness) ?
C e.g. pion-nucleon case
C
else if(is.eq.0.and.iabs(i2).le.maxbar)then
c... (anti-)N*,D*
call anndex(0,m1,i1,iz1,m2,i2,iz2,sqrts,sig,gam,
. maxbra,minnuc+1,maxdel,brtype,branres)
C
C Check if second particle is baryon
C (with amount one of netstrangeness) ?
C
else if(is.eq.1.and.iabs(i2).le.maxbar)then
c... (anti-)Y*
call anndex(0,m1,i1,iz1,m2,i2,iz2,sqrts,sig,gam,
. maxbrs1,minlam+1,maxsig,bs1type,branbs1)
C
C Check if second particle is baryon
C (with amount two of netstrangeness) ?
C
else if(is.eq.2.and.iabs(i2).le.maxbar)then
c... (anti-)X*
call anndex(0,m1,i1,iz1,m2,i2,iz2,sqrts,sig,gam,
. maxbrs2,mincas+1,maxcas,bs2type,branbs2)
C
C
else
sig=0d0
return
end if
C
C *************************************************************
Css Case 2 : one ingoing particle (resonance,..) --> 2-4 outgoing
Css particles (decay)
C
else ! decay !!!!!
i2=0
iz2=0
m2=0.d0
is=iabs(strit(i1))
c
if(iabs(i1).ge.minmes)then ! meson dec.
call anndex(1,m1,i1,iz1,m2,i2,iz2,sqrts,sig,gam,
. maxbrm ,minmes+1,maxmes,bmtype,branmes)
else if(is.eq.0)then ! n*,d,d*
call anndex(1,m1,i1,iz1,m2,i2,iz2,sqrts,sig,gam,
. maxbra,minnuc+1,maxdel,brtype,branres)
else if(is.eq.1)then !
call anndex(1,m1,i1,iz1,m2,i2,iz2,sqrts,sig,gam,
. maxbrs1,minlam+1,maxsig,bs1type,branbs1)
else if(is.eq.2)then
call anndex(1,m1,i1,iz1,m2,i2,iz2,sqrts,sig,gam,
. maxbrs2,mincas+1,maxcas,bs2type,branbs2)
else
write(6,*)'make22(anndex): s=',is,'not included'
stop
end if
C
C End of Cases 1 and 2 : annihilation/decay
C
end if
C ************************************************************
C
return
end
C####C##1#########2#########3#########4#########5#########6#########7##
subroutine anndex(io,m1,i1,iz1,m2,i2,iz2,sqrts,sig,gam,
& maxbr,mini,maxi,btype,branch)
c
cinput io : 0: annihilation; 1: decay
cinput m1 : mass of scattering/decaying particle 1
cinput i1 : ID of scattering/decaying particle 1
cinput iz1 : $2\cdot I_3$ of scattering/decaying particle 1
cinput m2 : mass of scattering particle 2
cinput i2 : ID of scattering particle 2
cinput iz2 : $2\cdot I_3$ of scattering particle 2
cinput sqrts : $sqrts{s}$ of collision; resonance mass for decay
coutput sig : resonance scattering cross section
coutput gam : width of the produced resonance
cinput maxbr : number of decay channels for particle class
cinput mini : smallest {\tt ityp} of particle class
cinput maxi : largest {\tt ityp} of particle class
cinput btype : array with exit channel definitions
cinput branch : array with branching ratios for final state
c
c
c {\tt anndex} performs meson-baryon and meson-meson annihilations
c as well as all meson and baryon resonance decays. In case of
c annihilations it returs the total resonance production cross section
c and the decay width of the resonance chosen as final state.
c The final state itself for both cases, annihilation and decay
c is returned via the {\tt newpart} common block. In the case
c of a decay the final state may consist of up to 4 particles.
c
c In {\tt anndex} the actual summation over Breit-Wigner formulas
c in the case of annihilations is performed.
c
c For decays the branch is choosen according to the mass dependent
c part of the decay width (call to {\tt fbrancx}); then the final
c state (which consists of particle {\tt ityp}, $2\cdot I_3$ and
c mass is generated and transferred to teh {\tt newpart} common
c block.
c
Ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
implicit none
include 'comres.f'
include 'comwid.f'
include 'newpart.f'
include 'options.f'
real*8 pi,cc,sqrts
parameter(pi=3.1415927,cc=0.38937966)
integer maxbr,mini,maxi,btype(4,0:maxbr)
real*8 branch(0:maxbr,mini:maxi)
integer icnt,is
integer io,i,j,i1,i2,iz1,iz2,itag,ii1,ii2
integer itn1,itnz1
real*8 m1,m2,prob(0:100),sum,sig,gam,cgk2
real*8 sigi(minnuc:maxmes),mmax,mmin,br,mmi1,mmi2,ppcm,gt
real*8 m,g,mo
c functions
real*8 fbrwig,pcms,massit
real*8 mminit,widit,fbrancx,ranf,fcgk,fwidth
real*8 fprwdt
integer jit,isoit,strit
if(io.eq.1)then
C
C
C one ingoing particle --> two,three,four outgoing particles
C
c... decays
do 3 i=0,maxbr
if(isoit(btype(1,i))+isoit(btype(2,i))+isoit(btype(3,i))+
& isoit(btype(4,i)).lt.iabs(iz1).or.
& m1.lt.mminit(btype(1,i))+mminit(btype(2,i))
& +mminit(btype(3,i))+mminit(btype(4,i)) )then
prob(i)=0.d0
else
prob(i)=fbrancx(i,iabs(i1),iz1,m1,branch(i,iabs(i1)),
& btype(1,i),btype(2,i),btype(3,i),btype(4,i))
endif
3 continue
icnt=0
c... find out branch = i
call getbran(prob,0,100,sum,0,maxbr,i)
ctp060202 9 call getbran(prob,0,100,sum,0,maxbr,i)
if(i.gt.maxbr)then
write(6,*)'anndex(dec): no final state found for:',i1,m1,iz1
write(6,*)'please check minimal masses: m1,m1min,m2min'
write(6,*)'and iso3 of decaying particle'
write(6,*)(prob(j),j=0,maxbr)
stop
end if
c... get itypes and set number of outgoing particles, prepare final state
nexit=2
itypnew(1)=btype(1,i)
itot(1)=isoit(itypnew(1))
itypnew(2)=btype(2,i)
itot(2)=isoit(itypnew(2))
itypnew(3)=btype(3,i)
if(itypnew(3).ne.0) then
itot(3)=isoit(itypnew(3))
pnew(5,3)=massit(itypnew(3))
c sidnew is used to set the lstcoll array
sidnew(3)=strcount
sidnew(2)=strcount ! correct here, set only for nexit > 2
nexit=nexit+1
endif
itypnew(4)=btype(4,i)
if(itypnew(4).ne.0) then
itot(4)=isoit(itypnew(4))
pnew(5,4)=massit(itypnew(4))
sidnew(4)=strcount
nexit=nexit+1
endif
if(nexit.gt.2) strcount=strcount+1
c check for some special cases involving decay of antibaryons and
c strange mesons
if(iabs(i1).ge.minmes.and.strit(i1).ne.0) then
do 41 j=1,nexit
if(strit(itypnew(j)).ne.0)then
c for anti-K* decays(mesons with one s-quark)
itypnew(j)=isign(itypnew(j),i1)
end if
41 continue
elseif(iabs(i1).lt.minmes) then
c the (anti-)baryon MUST always be the first outgoing particle
c -> conserve baryon-charge
itypnew(1)=isign(itypnew(1),i1)
do 42 j=2,nexit
if(strit(itypnew(j)).ne.0.and.i1.lt.0) then
itypnew(j)=(-1)*itypnew(j)
endif
42 continue
endif
c... get isopin-3 components
itag=-50
call isonew4(isoit(i1),iz1,itot,i3new,itag)
c write(6,*)'anndem:',iz3,iz1,iz2,'#',is3,is1,is2
c... get masses
if(widit(itypnew(1)).ge.1.d-4.and.
& widit(itypnew(2)).le.1.d-4)then
c... i1 is a broad meson
pnew(5,2)=massit(itypnew(2))
mmin=mminit(itypnew(1))
mo = pnew(5,2)
if(nexit.gt.2) then
do 39 j=3,nexit
mo=mo+pnew(5,j)
39 continue
endif
mmax=sqrts-mo
call getmas(massit(itypnew(1)),widit(itypnew(1)),itypnew(1)
& ,isoit(itypnew(1)),mmin,mmax,mo,pnew(5,1))
elseif(widit(itypnew(2)).ge.1.d-4
& .and.widit(itypnew(1)).le.1.d-4)then
c... i2 is a broad meson
pnew(5,1)=massit(itypnew(1))
mmin=mminit(itypnew(2))
mo = pnew(5,1)
if(nexit.gt.2) then
do 49 j=3,nexit
mo=mo+pnew(5,j)
49 continue
endif
mmax=sqrts-mo
call getmas(massit(itypnew(2)),widit(itypnew(2)),itypnew(2)
& ,isoit(itypnew(2)),mmin,mmax,mo,pnew(5,2))
elseif(widit(itypnew(1)).ge.1.d-4
& .and.widit(itypnew(2)).ge.1.d-4)then
c... i1&i2 are both broad
if(ranf(0).gt.0.5)then
mmin=mminit(itypnew(1))
mo=mminit(itypnew(2))
if(nexit.gt.2) then
do 59 j=3,nexit
mo=mo+pnew(5,j)
59 continue
endif
mmax=sqrts-mo
call getmas(massit(itypnew(1)),widit(itypnew(1)),
& itypnew(1),isoit(itypnew(1)),mmin,mmax,mo,pnew(5,1))
mmin=mminit(itypnew(2))
mo=pnew(5,1)
if(nexit.gt.2) then
do 69 j=3,nexit
mo=mo+pnew(5,j)
69 continue
endif
mmax=sqrts-mo
call getmas(massit(itypnew(2)),widit(itypnew(2)),
& itypnew(2),isoit(itypnew(2)),mmin,mmax,mo,pnew(5,2))
else ! of ranf.gt.0.5
mmin=mminit(itypnew(2))
mo=mminit(itypnew(1))
if(nexit.gt.2) then
do 79 j=3,nexit
mo=mo+pnew(5,j)
79 continue
endif
mmax=sqrts-mo
call getmas(massit(itypnew(2)),widit(itypnew(2)),
& itypnew(2),isoit(itypnew(2)),mmin,mmax,mo,pnew(5,2))
mmin=mminit(itypnew(1))
mo=pnew(5,2)
if(nexit.gt.2) then
do 89 j=3,nexit
mo=mo+pnew(5,j)
89 continue
endif
mmax=sqrts-mo
call getmas(massit(itypnew(1)),widit(itypnew(1)),
& itypnew(1),isoit(itypnew(1)),mmin,mmax,mo,pnew(5,1))
endif
c none are broad
else
pnew(5,2)=massit(itypnew(2))
pnew(5,1)=massit(itypnew(1))
end if
mmax=0.d0
do 99 j=1,nexit
mmax=mmax+pnew(5,j)
99 continue
if(sqrts.le.mmax)then
write(6,*)' *** error(anndex): treshold violated',sqrts-mmax
stop
end if
C
C
C two ingoing particles --> one outgoing particle (resonance)
C
C i.e. (i0=0)
else
c.... collisions: find in-branch = j
sig=0.0
gam=0.0
C
ii1=i1
ii2=i2
c for strange - nonstrange meson-meson scattering: strip sign
is=iabs(strit(i1)+strit(i2))
if(is.ne.0.and.iabs(i1).ge.minmes.and.iabs(i2).ge.minmes) then
ii1=iabs(i1)
ii2=iabs(i2)
endif
c for meson baryon, strip sign of baryon
if(iabs(i2).le.maxbar) then
ii2=iabs(i2)
endif
c
call getobr(btype,0,maxbr,ii1,ii2,j)
if(j.eq.-99)return
mmi1=mminit(i1)
mmi2=mminit(i2)
c
C next line outside the loop (compare post-QM-version: inside the loop)
C
C
ppcm=pcms(sqrts,m1,m2)
C
C Loop over different branches (resonances...)
C
do 88 i=mini,maxi
sigi(i)=0.d0
br=branch(j,i)
gt=widit(i)
if(br*gt.lt.1d-4)goto 88
cgk2=fcgk(i1,i2,iz1,iz2,i)
if(br*cgk2.gt.0d0.and.sqrts.gt.mmi1+mmi2+1d-2.and.
& ppcm.gt.1d-2)then
C
C
br=fprwdt(j,i,iz1+iz2,sqrts)/fwidth(i,iz1+iz2,sqrts)
m=dabs(sqrts)
g=fwidth(i,iz1+iz2,m)
sigi(i)=dble(jit(i)+1)
/ /dble((jit(i1)+1)*(jit(i2)+1))
* *pi/ppcm**2*br
* *g*g/((m-massit(i))**2+g*g/4d0)*cgk2*cc
end if
if(sigi(i).gt.1e10)then
write(6,*)' ***error(anndec) cross section too high '
write(6,*)'anndex(ann):',i,
, br,cgk2,fbrwig(i,iz1+iz2,sqrts,1),
, 1/pcms(sqrts,m1,m2),sigi(i)
write(6,*)m1,m2,sqrts
write(6,*)i1,i2,i
write(6,*)iz1,iz2,iz1+iz2
end if
C
C
88 continue
c... find outgoing resonance
call getbran(sigi,minnuc,maxmes,sig,mini,maxi,itn1)
ctp060202 108 call getbran(sigi,minnuc,maxmes,sig,mini,maxi,itn1)
if(sig.ge.1d-10)then
itnz1=iz1+iz2
gam=fwidth(itn1,itnz1,sqrts)
end if
c copy created resonance into newpart arrays
itypnew(1)=itn1
i3new(1)=itnz1
pnew(5,1)=sqrts
if(iabs(i1).ge.minmes.and.iabs(i2).ge.minmes) then
if(iabs(strit(i1)+strit(i2)).ne.0) then
itypnew(1)=isign(itypnew(1),i1*i2)
endif
else
itypnew(1)=isign(itypnew(1),i2)
endif
ctp060202 3333 continue
C
C End of two Cases (annihilation/decay)
C
C
end if !dec/ann
return
end
C####C##1#########2#########3#########4#########5#########6#########7##
real*8 function fbrwig(i,iz,mi,bit)
c
cinput i : resonance ID
cinput iz : $2\cdot I_3$ of resonance
cinput mi : mass of resonance
cinput bit : sign is used as option to toggle between fixed and m.dep. widths
c
c {\tt fbrwig} returns a normalized Breit-Wigner Function.
c Note, that the normalization actually only holds true for
c fixed decay widths. {\tt fbrwig}, however, uses per default
c a mass dependent width. You should divide by
c {\tt bwnorm} to obtain normalized Breit-Wigners for mass dependent
c widths also. For {\tt bit} < 0 a fixed width is used.
c
cccccCcc1ccccccccc2ccccccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc
implicit none
integer i,iz,bit
real*8 pi,mi,g,fwidth,massit,widit
real*8 f,e,m0,g1,g2
parameter(pi=3.1415927)
f(e,m0,g1,g2)=0.5/pi*g1/((e-m0)**2+0.25*g2**2)
if(bit.lt.0)then
g=widit(i)
fbrwig=f(mi,massit(i),g,g)
else
g=fwidth(i,iz,mi)
fbrwig=f(mi,massit(i),g,g)
end if
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine getbran(x,dmin,dmax,sumx,nmin,nmax,i)
c
c
cinput x : vector containing weights, dimension is {\tt x(dmin:dmax)}
cinput dmin : lower dimension of {\tt x}
cinput dmax : upper dimension of {\tt x}
coutput sumx : sum of elements of {\tt x} from {\tt nmin} to {\tt nmax}
cinput nmin : lower boundary for {\tt getbran} operation
cinput nmax : upper boundary for {\tt getbran} operation
coutput i : index of element which has been choosen randomly
c
c {\tt getbran} takes a vector of weights or probabilities
c {\tt x(dmin:dmax)} and sums up the elements from
c {\tt nmin} to {\tt nmax}. It then chooses randomly an element {\tt i}
c between {\tt nmin} and {\tt nmax}. The probability of
c choosing {\tt i} depends on the weights contained in {\tt x}.
c
c \danger{ {\tt i} will be undefined if {\tt sum} is less or
c equal to zero}
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
implicit none
integer i,j,nmin,nmax,dmin,dmax
real*8 x(dmin:dmax),sumx,ranf,rx,cut
parameter (cut=1d-20)
sumx=0D0
do 10 j=nmin,nmax
sumx=sumx+x(j)
10 continue
if (sumx.lt.cut) then
i=nmax+1
return
endif
rx=sumx*ranf(0)
do 20 j=nmin,nmax
if (rx.le.x(j)) then
i=j
return
endif
rx=rx-x(j)
20 continue
if (abs(rx).lt.1D-10) then
i=nmax
return
endif
call error ('getbran','no channel found',rx,3)
end
C####C##1#########2#########3#########4#########5#########6#########7##
subroutine getobr(x,dmin,dmax,i1,i2,i)
c
c
cinput x : array, either {\tt brtype, bmtype, bs1type} or {\tt bs2type}
cinput dmin : lower dimension of {\tt x(4,dmin:dmax)}
cinput dmin : upper dimension of {\tt x(4,dmin:dmax)}
cinput i1 : ID of first incoming particle
cinput i2 : ID of second incoming particle
coutput i : index of decay branch into {\tt i1} and {\tt i2}
c
c {\tt getobr} returns the index of the decay branch for the
c exit channel $B^* \rightarrow$ {\tt i1} + {\tt i2}
c from one of the arrays
c {\tt brtype, bmtype, bs1type} or {\tt bs2type}. This index
c is needed for the calculation of the cross section
c {\tt i1} + {\tt i2} $\rightarrow B^*$.
c
cccccCcc1ccccccccc2ccccccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc
implicit none
integer i,j,i1,i2,dmin,dmax,x(4,dmin:dmax)
do 108 j=dmin,dmax
if((x(1,j).eq.i1.and.x(2,j).eq.i2.and.x(3,j).eq.0).OR.
& (x(1,j).eq.i2.and.x(2,j).eq.i1.and.x(3,j).eq.0))then
i=j
return
end if
108 continue
i=-99
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine normit (sigma,isigline)
c
c Revision : 1.0
c
cinput sigma : vector with all (partial) cross sections
coutput sigma : unitarized vector with cross sections
cinput isigline : process class of cross sections
c
c {\tt normit} unitarizes the cross sections contained in the
c {\tt sigma} array. The total cross section is stored in
c {\tt sigma(0)}. The partial cross sections are unitarized
c (rescaled) such, that their sum adds up to the total cross
c section. Confidence levels can be assigned to different
c partial cross sections indicating whether they may be rescaled
c (i.e. if they are not well known) or whether they must not
c be rescaled (i.e. because they have been fitted to experimental
c data).
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
implicit none
include 'options.f'
include 'comres.f'
real*8 sigma(0:maxpsig)
integer isigline
integer i, npsig, restart
integer uncert(1:maxpsig)
real*8 diff, sumpart, gsumpart
real*8 newsig(0:maxpsig)
c get the number of channels
npsig=sigmaln(1,1,isigline)
c normalize only if sigtot is not given by the sum of sigpart
if (sigmaln(2,1,isigline).gt.0) then
c copy array
do 10 i=1,npsig
uncert(i)=sigmaln(i+2,2,isigline)
10 continue
100 restart=0
sumpart=0
gsumpart=0
c calculate the sum of all sigpart
do 20 i=1,npsig
sumpart=sumpart+sigma(i)
gsumpart=gsumpart+sigma(i)*uncert(i)
20 continue
c difference between sigtot and the sum of sigpart
diff=sigma(0)-sumpart
c if all channels are exactly zero, there must be an error in blockres.f!
if (sumpart.eq.0.0) then
write (6,*) 'normit: Error! sumpart.eq.0'
c stop
return
endif
if (gsumpart.eq.0.0) then
do 50 i=1,npsig
c now all channels can be modified
if (uncert(i).eq.0) then
c write (6,*) 'modify channel',i
uncert(i)=1
endif
50 continue
c restart calculation
goto 100
endif
do 60 i=1,npsig
c rescale channels
newsig(i)=sigma(i)+uncert(i)*diff*sigma(i)/gsumpart
c if a channel is negative...
if (newsig(i).lt.0) then
c set it to zero and restart
sigma(i)=0.0
restart=1
endif
60 continue
if (restart.eq.1) goto 100
c copy new values to sigma
do 70 i=1,npsig
sigma(i)=newsig(i)
70 continue
endif
if (CTOption(7).eq.1.and.sigma(2).gt.1d-10) then
sigma(1)=0d0
endif
if (CTOption(7).eq.-1) then
do 80 i=2,npsig
sigma(i)=0d0
80 continue
end if
return
end
C####C##1#########2#########3#########4#########5#########6#########7##
real*8 function fwidth(ir,izr,m)
c
cinput ir : resonance ID
cinput izr : $2\cdot I_3$ of resonance
cinput m : mass of resonance
c
c {\tt fwidth} returns the mass-dependent total decay width
c of the resonance {\tt ir}.
c
c
C####C##1#########2#########3#########4#########5#########6#########7##
implicit none
include 'comres.f'
include 'comwid.f'
include 'options.f'
integer i,ir,izr,mm,mp,ires
real*8 gtot,m,widit,splint
real*8 minwid, fprwdt
if (CTOption(1).ne.0) then
fwidth=widit(ir)
return
endif
if (wtabflg.gt.0.and.CTOption(33).eq.0) then
ires=iabs(ir)
minwid=min(widit(ir),1D-8)
if (ires.ge.minbar.and.ires.le.maxbar) then !baryons
c widths are continued horicontally outside the spline region
if(m.le.maxtab2)then
fwidth=max(splint(tabx(1),fbtaby(1,ires,1),
. fbtaby(1,ires,2),widnsp,m),minwid)
else
fwidth=max(splint(tabx(1),fbtaby(1,ires,1),
. fbtaby(1,ires,2),widnsp,maxtab2),minwid)
endif
else if (ires.ge.minmes.and.ires.le.maxmes) then !mesons
c widths are continued horicontally outside the spline region
if(m.le.maxtab2)then
fwidth=max(splint(tabx(1),fmtaby(1,ires,1),
. fmtaby(1,ires,2),widnsp,m),minwid)
else
fwidth=max(splint(tabx(1),fmtaby(1,ires,1),
. fmtaby(1,ires,2),widnsp,maxtab2),minwid)
endif
else
write (6,*) '*** error(fwidth) wrong itype:',ir
fwidth=0
endif
else
call brange(ir,mm,mp)
gtot=0d0
if(mp.gt.0)then
do 27 i=mm,mp
gtot=gtot+fprwdt(i,ir,izr,m)
27 continue
end if
fwidth=gtot !*widit(ir)
end if
return
end
C####C##1#########2#########3#########4#########5#########6#########7##
real*8 function fprwdt(i,ir,izr,mi)
c
cinput i : decay branch
cinput ir : resonance ID
cinput izr : $2\cdot I_3$ of resonance
cinput mi : mass of resonance
c
c {\tt fprwdt} returns the mass dependent partial decay width
c of the decay channel {\tt i}.
c
C####C##1#########2#########3#########4#########5#########6#########7##
implicit none
real*8 m,br,bi,mir,g,mi
integer i,ir,izr,i1,i2,i3,i4
real*8 widit,fbrancx,mminit,massit
call b3type(ir,i,bi,i1,i2,i3,i4)
m=dabs(mi)
g=0d0
mir=massit(ir)
if(bi.gt.1d-9.and.mir.gt.mminit(i1)+mminit(i2))then
br=fbrancx(i,ir,izr,m,bi,i1,i2,i3,i4)
c write(6,*)' ',bi,gi,widit(ir)
g=br*widit(ir)
end if
fprwdt=g
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
real*8 function fbrancx(i,ir,izr,em,bi,b1,b2,b3,b4)
c
cinput i : decay branch
cinput ir : ID of resonance
cinput em : actual mass of resonance
cinput bi : branching ration at peak
cinput b1 : itype of 1st outgoing particle
cinput b2 : itype of 2nd outgoing particle
cinput b3 : itype of 3rd outgoing particle
cinput b4 : itype of 4th outgoing particle
c
c {\tt fbrancx} returns the mass dependent branching ratio for
c the decay channel {\tt i} of resonance {\tt ir}. This
c branching ratio is NOT normalized. To extract the mass dependent
c decay width, use {\tt fprwdt}.
c
c {\tt fbrancx} =$
c \left( \Gamma^{i,j}_{R} \frac{M_{R}}{M}
c \left( \frac{\langle p_{i,j}(M) \rangle}
c {\langle p_{i,j}(M_{R}) \rangle} \right)^{2l+1}
c \frac{1.2}{1+ 0.2
c \left( \frac{\langle p_{i,j}(M) \rangle}
c {\langle p_{i,j}(M_{R}) \rangle} \right)^{2l} }
c \right) $
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
implicit none
real*8 kdiv1,kdiv2,em,b,mmin,mn,m1m,m2m
real*8 bi,minwid
real*8 fbran,splint,splintth,pmean
integer i,ires,ir,izr,b1,b2,b3,b4
include 'comres.f'
include 'comwid.f'
include 'options.f'
real*8 mminit,massit
integer isoit,flbr,ipwr,ipwr1
ires=iabs(ir)
if(iabs(izr).gt.isoit(ires))then
fbrancx=0d0
return
end if
if(CTOption(8).ne.0)then
fbrancx=bi
return
end if
m1m=mminit(b1)
m2m=mminit(b2)
c in case of three or four particle decays put masses in m2m
if(b3.ne.0) m2m=m2m+mminit(b3)
if(b4.ne.0) m2m=m2m+mminit(b4)
mn=massit(ires) ! nominal mass
mmin= m1m+m2m ! minimal mass of resonance
if (wtabflg.ge.2.and.CTOption(33).eq.0) then
minwid=min(fbran(i,ires),1D-8)
if (ires.ge.minbar.and.ires.le.maxbar) then !baryons
c branching ratios are continued horicontally outside the spline region
if(em.le.maxtab2)then
b=max(splintth(tabx,pbtaby(1,1,ires,i),
. pbtaby(1,2,ires,i),widnsp,em,mmin),minwid)
else
b=max(splintth(tabx,pbtaby(1,1,ires,i),
. pbtaby(1,2,ires,i),widnsp,maxtab2,mmin),minwid)
endif
else if (ires.ge.minmes.and.ires.le.maxmes) then !mesons
if (em.le.maxtab2) then
c branching ratios are continued horicontally outside the spline region
b=max(splint(tabx,pmtaby(1,1,ires,i),
. pmtaby(1,2,ires,i),widnsp,em),minwid)
else
b=max(splint(tabx,pmtaby(1,1,ires,i),
. pmtaby(1,2,ires,i),widnsp,maxtab2),minwid)
endif
else
write (6,*) '*** error(fbrancx) wrong id:',ir
b=0
endif
else
b=0d0
if (bi.gt.0.and.em.gt.mmin.and.mn.gt.mmin) then
ipwr=flbr(i,ires)
ipwr1=ipwr+1
c determine expectation values of outgoing masses
c call of pmean with -99 instead of iso3 to ensure usage of fixed
c resonance widths: 5% error, but avoids recursion via call
c to fwidth from pmean
if(CTOption(33).ne.0)then
kdiv1=pmean(em,b1,-99,b2,-99,b3,-99,b4,-99,ipwr1)/
& pmean(mn,b1,-99,b2,-99,b3,-99,b4,-99,ipwr1)
kdiv2=pmean(em,b1,-99,b2,-99,b3,-99,b4,-99,ipwr)/
& pmean(mn,b1,-99,b2,-99,b3,-99,b4,-99,ipwr)
else
kdiv1=pmean(em,b1,isoit(b1),b2,isoit(b2),
& b3,isoit(b3),b4,isoit(b4),ipwr1)/
& pmean(mn,b1,isoit(b1),b2,isoit(b2),
& b3,isoit(b3),b4,isoit(b4),ipwr1)
kdiv2=pmean(em,b1,isoit(b1),b2,isoit(b2),
& b3,isoit(b3),b4,isoit(b4),ipwr)/
& pmean(mn,b1,isoit(b1),b2,isoit(b2),
& b3,isoit(b3),b4,isoit(b4),ipwr)
end if
b=bi*mn/em*kdiv1*1.2/(1.+0.2*kdiv2)
else
b=0.
end if
end if
fbrancx=b
return
end