-
ralfulrich authored
rename framweork/sequence int framework/process remove all dynamic build files from main corsika directory, move into src re-added main/shower fixed wrong file locations Update corsika.hpp. Re-defining the version macros. Update corsika.hpp Update CONTRIBUTING.md Update MAINTAINERS.md Deleted FIXME.md Deleted AUTHORS Update COLLABORATION_AGREEMENT.md Update .clang-format Update .gitlab-ci.yml added packages: proposal, conex, cnpy, spdlog, lcov, corsika-data prevent in-source builds fotran language flags, sections build types fixed stack_example rename directories corsika8 interface target fixed few compilation, dependency issues in new system clang error simpler cmake files, using functions, re-introduce run_examples target clang format copyright notices fixed CI script, clang-format simplify ctest output better interfaces for urqmd and qgsjII, conex updated conex, data cmake integration build system cmake updates also updated corsika.hpp added earth radius moved static_pow from corsika::units::si::detail to corsika::units updated units::si namespaces throughout the project No python jobs removed a lot of stray using namespaces and corsika::units::si in entire codebase [refactory-2020] stack implementations: updated [refactory-2020] stack implementations: updated [refactory-2020] stack implementations [refactory-2020] stack implementations: !290 headers [refactory-2020] stack implementations: !290 headers
ralfulrich authoredrename framweork/sequence int framework/process remove all dynamic build files from main corsika directory, move into src re-added main/shower fixed wrong file locations Update corsika.hpp. Re-defining the version macros. Update corsika.hpp Update CONTRIBUTING.md Update MAINTAINERS.md Deleted FIXME.md Deleted AUTHORS Update COLLABORATION_AGREEMENT.md Update .clang-format Update .gitlab-ci.yml added packages: proposal, conex, cnpy, spdlog, lcov, corsika-data prevent in-source builds fotran language flags, sections build types fixed stack_example rename directories corsika8 interface target fixed few compilation, dependency issues in new system clang error simpler cmake files, using functions, re-introduce run_examples target clang format copyright notices fixed CI script, clang-format simplify ctest output better interfaces for urqmd and qgsjII, conex updated conex, data cmake integration build system cmake updates also updated corsika.hpp added earth radius moved static_pow from corsika::units::si::detail to corsika::units updated units::si namespaces throughout the project No python jobs removed a lot of stray using namespaces and corsika::units::si in entire codebase [refactory-2020] stack implementations: updated [refactory-2020] stack implementations: updated [refactory-2020] stack implementations [refactory-2020] stack implementations: !290 headers [refactory-2020] stack implementations: !290 headers
string.f 57.63 KiB
c $Id: string.f,v 1.18 2003/05/02 11:06:56 weber Exp $
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine stringdec(ityp,iz2,smass,part,ident2,npart)
c
cinput smass : Stringmass
cinput ityp : Particle ID
cinput iz2 : Isospin$_3\cdot 2$
c
coutput part : 4-momenta, 4-position, masses (array)
coutput ident2 : ityp, iz2 (array)
coutput npart : number of outgoing particles
c
c This subroutine performs string fragmentation.
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
implicit real*8 (a-h,o-z)
implicit integer (i-n)
PARAMETER(MXPTCL=200)
COMMON/PARTCL/ PPTCL(9,MXPTCL),nptcl,IDENT(MXPTCL),IDCAY(MXPTCL)
include 'comstr.f'
dimension part(9,mxptcl)
dimension ident2(2,mxptcl)
c we call the translation routine
call ityp2id(ityp,iz2,ifa,ifb)
goto 1
cspl... string fragmentation called with quark id's and energy as arguments
entry qstring(ifanew,ifbnew,smass,part,ident2,npart)
ifa=ifanew
ifb=ifbnew
1 continue
smem=smass
c here we call the fragmentation routine. the produced hadrons and their
c properties are returned via the pptcl- and ident-array in the common-block
call string(ifa,ifb,smass)
c here the array pptcl has been filled with nptcl entries (1->nptcl)
c now we translate to uqmd and shift the pptcl- and ident-info to the
c corresponding part- and ident2-arrays of the newpart-common-block:
npart=nptcl
do 2 i=1,nptcl
call id2ityp(ident(i),pptcl(5,i),itypout,iz2out)
ident2(1,i)=itypout
ident2(2,i)=iz2out
smem=smem-pptcl(4,i)
do 3 j=1,9
part(j,i)=pptcl(j,i)
3 continue
2 continue
c check for energy conservation:
ctp060926 if(abs(smem).gt.1.0D-5)then
ctp060926 write(*,*)'! stringdec: energy difference=',smem
ctp060926 write(*,*)'ifa,ifb,smass=',ifa,ifb,smass
ctp060926 endif
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine strini
c
c output : via common blocks
c
c {\tt strini} calculates mixing angles for the meson-multiplettsc
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
implicit real*8 (a-h,o-z)
implicit integer (i-n)
include 'options.f'
include 'comres.f'
include 'comstr.f'
real*8 m3
real*8 massit
integer jit
c mixing angles of meson multiplets according to flavor SU(3) quark model:
cspl-0795 these parameters assign the pure u/ubar,d/dbar,s/sbar states
c (e.g. 110,220,330) to the physical particles according to the su(3)
c quark model. The flavor mixing angles are chosen according to quadratic
c Gell-Mann-Okubo mass formula (Review of Particle Properties,
c Phys. Rev D50 (1994) 1319). For the scalar mesons this formula is not
c applicable. We assume an ideal mixing angle (tan(theta)=1/sqrt(2)).
c
c pseudoscalar: theta=-10 deg
c vector : theta= 39 deg
c pseudovector: theta= 51 deg
c tensor : theta= 28 deg
c
real*8 mixang(njspin)
c ideal mixing angles assumed for the last three multiplets
data mixang/-10d0,39d0,35.3d0,51d0,28d0,35.3d0,35.3d0,35.3d0/
c
pi = 4d0*datan2(1d0,1d0)
c calculate 'singlet shift probabilities', e.g. a 11 (u-ubar) state
c can be changed to a 22 (d-dbar) or 33 (s-sbar) state with a certain
c probability. THEN they can be identified with physical hadrons!
do 3 i=1,njspin
mixang(i)=mixang(i)/36d1*2d0*3.1416d0
PMIX1S(1,i)=(dcos(mixang(i))/sqrt(6d0)
& -dsin(mixang(i))/sqrt(3d0))**2
PMIX1S(2,i)=PMIX1S(1,i)
PMIX1S(3,i)=(-(dcos(mixang(i))*2d0/dsqrt(6d0))
& -dsin(mixang(i))/dsqrt(3d0))**2
PMIX2S(1,i)=0.5d0
PMIX2S(2,i)=0.5d0
PMIX2S(3,i)=1d0
ce calculate probabilities of the meson multipletts
ce according to parm=(spin degeneracy)/(average mass) *ctp(50 ff.)
parm(i)=0d0
m3=0d0
do 102 j=0,3
itp=mlt2it(4*(i-1)+j+1)
m3=m3+massit(itp)
jpc=jit(itp)/2
102 continue
parm(i)=parm(i)+(2*jpc+1)/m3*4*CTParam(49+i)
c the mixing-angles are the same for 'string' and 'cluster':
do 2 k=1,3
PMIX1C(k,i)=PMIX1S(k,i)
PMIX2C(k,i)=PMIX2S(k,i)
c write(6,*)'#',pmix1c(k,i)
2 continue
3 continue
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
SUBROUTINE GAUSPT(PT0,SIGQT)
ccinput sigqt : Width of Gaussian
c
coutput pt0 : transverse momentum
c
C generate pt with Gaussian
c distribution $\propto pt \exp(-pt^2/sigqt^2)$
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
implicit none
real*8 pt0,sigqt,rnd,ranf
RND=ranf(0)
PT0=SIGQT*SQRT(-DLOG(1.d0-RND))
RETURN
END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
SUBROUTINE FLAVOR(ID,IFL1,IFL2,IFL3,JSPIN)
c
cinput ID : quarkcode
c
coutput ifl1 : single quarks id
coutput ifl2 : single quarks id
coutput ifl3 : single quarks id
coutput jspin : spin id
c
C THIS SUBROUTINE UNPACKS THE IDENT CODE ID=+/-IJKL
c
C -MESONS:
C I=0, J<=K, +/- IS SIGN FOR J,
C ID=110 FOR PI0, ID=220 FOR ETA, ETC.
c
C -BARYONS:
C I<=J<=K IN GENERAL,
C J<I<K FOR SECOND STATE ANTISYMMETRIC IN (I,J), EG. L = 2130
c
C -DIQUARKS:
C ID=+/-IJ00, I<J FOR DIQUARK COMPOSED OF I,J.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
implicit real*8 (a-h,o-z)
implicit integer (i-n)
IDABS=IABS(ID)
c extract the single (anti-)quark ids from the hadron id:
I=IDABS/1000
J=MOD(IDABS/100,10)
K=MOD(IDABS/10,10)
JSPIN=MOD(IDABS,10)
c diquarks:
IF(ID.NE.0.AND.MOD(ID,100).EQ.0) GO TO 300
c quarks oder so:
IF(J.EQ.0) GO TO 200
c mesonen:
IF(I.EQ.0) GO TO 100
C..BARYONS:
C..ONLY X,Y BARYONS ARE QQX, QQY, Q=U,D,S.
IFL1=ISIGN(I,ID)
IFL2=ISIGN(J,ID)
IFL3=ISIGN(K,ID)
RETURN
C MESONS
100 CONTINUE
IFL1=0
IFL2=ISIGN(J,ID)
IFL3=ISIGN(K,-ID)
RETURN
200 CONTINUE
IFL1=0 IFL2=0
IFL3=0
JSPIN=0
return
300 IFL1=ISIGN(I,ID)
IFL2=ISIGN(J,ID)
IFL3=0
JSPIN=0
RETURN
END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
SUBROUTINE STRING(IFL1,IFL2,AMSTR)
c
cinput amstr : stringmass
cinput ifl1 : leading (di)quark (along +Z)
cinput ifl2 : leading (di)quark
c
c output : produced particles via common block ({\tt pptcl})
c
C Hadron production via string fragmentation. masses acc. to Breit
c Wigner distr., incl. production of mesonic and baryonic resonances
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
implicit real*8 (a-h,o-z)
implicit integer (i-n)
COMMON/COMTRY/ NTRIES
PARAMETER(MXPTCL=200)
COMMON/PARTCL/ PPTCL(9,MXPTCL),nptcl,IDENT(MXPTCL),IDCAY(MXPTCL)
include 'comstr.f'
DIMENSION PX1L(2),PY1L(2),PX1(2),PY1(2),PMTS(2),W(2),IFL(2)
LOGICAL DIQBR,SPINT,BACK
DIMENSION VC(3)
DIMENSION PPTL(5,MXPTCL),PPTR(5,MXPTCL),
*IDENTL(MXPTCL),IDENTR(MXPTCL)
COMMON/KAPPA/ XAP
COMMON/COLRET/ LRET
LOGICAL LRET
COMMON/CONSTI/ CONSTI
LOGICAL CONSTI
LOGICAL leading
include 'options.f'
IFLL=0
PXL=0
PYL=0
c..strange quark and charm quark probabilities:
prbs=ctparam(6)
prbc=ctparam(7)
cspl the diquark-suppresion parameter is reduced for small
c string masses (finite size effect) see A. Jahns diploma thesis
if(amstr.le.2.8d0)then
pbars=0.d0
elseif (amstr.le.5d0)then
pbars=((amstr-2.8d0)/2.2d0)**3*ctparam(8)
else
pbars=ctparam(8)
endif
c..some parameters (see input.f)
dmas=ctparam(9)
dmass=ctparam(10)
pardbs=ctparam(25)
sigqts=ctparam(42)C
C PRBS STRANGENESS SUPPRESSION PARAMETER
C PRBC CHARM SUPPRESSION PARAMETER
PU=1./(2.+PRBS+PRBC)
C
LRET=.FALSE.
C
NREP = 0
DIQBR=.TRUE.
NFIX=0
BACK=.TRUE.
c.. here starts everything (if the string break up didn't work
c.. we start again at this point)
100 I=NFIX
ilead=0
NPR=0
NPL=0
NPTCL=NFIX
c.. nrep=number of tries to break string
NREP=NREP+1
IF(NREP.LT.NTRIES) GO TO 102
LRET=.TRUE.
ctp060926 write(*,*)'!! STRING: no fragmentation.'
return
102 CONTINUE
IFL(1)=IFL1
IFL(2)=IFL2
DO 110 J=1,2
110 W(J)=AMSTR
DO 120 J=1,2
PX1L(J)=0.
PY1L(J)=0.
PX1(J)=0.
120 PY1(J)=0.
C WILL THERE BE ONLY ONE BREAK OR NOT ?
SPINT=.TRUE.
KSPIN=1
IF(MOD(IFL(1),100).EQ.0.AND.MOD(IFL(2),100).EQ.0) GOTO 131
IDR=IDPARS(IFL(1),IFL(2),SPINT,KSPIN)
WEND=(AMASS(IDR)+DMAS)**2
GO TO 151
131 IFCN=1
IF(ranf(0).GT.0.5) IFCN=2
IFLC1=IFCN
IF(IFL(1).LT.0) IFLC1=-IFCN
IDR1=IDPARS(IFL(1),IFLC1,SPINT,KSPIN)
IDR2=IDPARS(IFL(2),-IFLC1,SPINT,KSPIN)
WEND=(AMASS(IDR1)+AMASS(IDR2)+DMAS)**2
151 SPINT=.FALSE.
KSPIN=0
c.. only one break goto 225 is the end of the fragmentation
IF(W(1)*W(2).LE.WEND) GO TO 225
c..the main iteration loop for the fragmentation:
130 I=I+1
IF(I.GT.MXPTCL) GO TO 9999
C CHOOSE SIDE OF BREAK
JSIDE=INT(1.+2.*ranf(0))
IF(JSIDE.EQ.1) NPR=NPR+1
IF(JSIDE.EQ.2) NPL=NPL+1
IF(NPR.GT.MXPTCL.OR.NPL.GT.MXPTCL) GO TO 9999
C IS IFL(JSIDE) A QUARK OR A DIQUARK ? (di-quark-->150)
IF(MOD(IFL(JSIDE),100).EQ.0) GO TO 150
C IFL(JSIDE) IS A QUARK
C NOW WE SELECT Q,QBAR PAIR OR QQ,QQBAR PAIR
DIQBR=.FALSE.
DRND=ranf(0)
c.. do a qq-qqbar pair with certain prob.
IF(DRND.LT.PBARS) GO TO 140C Q,QBAR PAIR
IFLN=ISIGN(IFLAV(PU,PRBS),-IFL(JSIDE))
GO TO 200
C QQ,QQBAR PAIR
140 IQ1=IFLAV(PU,PRBS)
IQ2=IFLAV(PU,PRBS)
c.. no single-strange diquarks (us,ds)!
cblubb if(max0(iq1,iq2).eq.3.and.min0(iq1,iq2).lt.3)goto 140
c.. suppr. double strange di-quarks(ss) with certain prob. (acc. to ctp 29)
if((IQ1.eq.3.and.IQ2.eq.3)
& .and.ranf(0).gt.CTParam(29))goto 140
IF(IQ1.LE.IQ2) GO TO 145
ISWAP=IQ1
IQ1=IQ2
IQ2=ISWAP
145 IFQQ=1000*IQ1+100*IQ2
IFLN=ISIGN(IFQQ,IFL(JSIDE))
GO TO 200
c..the di-quark part:
C IFL(JSIDE) IS A DIQUARK
C CAN DIQUARK BREAK OR NOT
150 DRND=ranf(0)
IF(DRND.LE. PARDBS) GO TO 190
C DIQUARK BREAK (prob. in PARDBS)
CALL FLAVOR(IFL(JSIDE),IFLD1,IFLD2,IFLD3,JSPIN)
IFLL=IFLD1
IFL(JSIDE)=IFLD2
DRND=ranf(0)
IF(DRND.GE.PARQLS) GO TO 160
IFLL=IFLD2
IFL(JSIDE)=IFLD1
160 DIQBR=.TRUE.
C LEADING QUARK TRANSVERSE MOMENTUM
CALL GAUSPT(PTL0,SIGQTS)
PHI=2.*PI*ranf(0)
PXL=PTL0*COS(PHI)
PYL=PTL0*SIN(PHI)
PX1L(JSIDE)=PX1(JSIDE)
PY1L(JSIDE)=PY1(JSIDE)
PX1(JSIDE)=-PXL
PY1(JSIDE)=-PYL
C Q,QBAR PAIR
IFLN=ISIGN(IFLAV(PU,PRBS),-IFL(JSIDE))
GO TO 200
C DIQUARK DOES NOT BREAK
C Q,QBAR PAIR
190 IFLN=ISIGN(IFLAV(PU,PRBS),IFL(JSIDE))
DIQBR=.FALSE.
C IDENT,MASS AND TRANSVERSE MOMENTUM OF PARTICLE
200 IDENT(I)=IDPARS(IFL(JSIDE),IFLN,SPINT,KSPIN)
PPTCL(5,I)=AMASS(IDENT(I))
SIGQTSN=SIGQTS
IF(MOD(IFLN,100).EQ.0) SIGQTSN=sigqts*ctparam(38)
if(iabs(ifln).eq.3.or.iabs(ifl(jside)).eq.3)
& sigqtsn=sigqts*ctparam(39)
CALL GAUSPT(PT2,SIGQTSN)
c no pt for leading hadron:
leading=.false.
if((JSIDE.EQ.1.and.NPR.eq.1).and.
& (abs(IDENT(I)).eq.1120 .or.abs(IDENT(I)).eq.1220) )then
c & (abs(IDENT(I)).ge.1110))then
leading=.true.
cblu pt2=pt2/2d0
ilead=ilead+1
endif
if((JSIDE.EQ.2.and.NPL.eq.1).and.
& (abs(IDENT(I)).eq.1120 .or.abs(IDENT(I)).eq.1220) )then
c & (abs(IDENT(I)).ge.1110))then leading=.true.
cblu pt2=pt2/2d0
ilead=ilead+1
endif
c..transverse momentum choosen for the newly produced hadron
PHI=2.*PI*ranf(0)
PX2=PT2*COS(PHI)
PY2=PT2*SIN(PHI)
PPTCL(1,I)=PX1(JSIDE)+PX2
PPTCL(2,I)=PY1(JSIDE)+PY2
C GENERATE Z-momentum
PMTS(3-JSIDE)=AMASS(IABS(IFL(3-JSIDE)))**2
PTS=PPTCL(1,I)**2+PPTCL(2,I)**2
PMTS(JSIDE)=PPTCL(5,I)**2+PTS
IF(PMTS(JSIDE)+PMTS(3-JSIDE).GE.PARRS*W(1)*W(2)) GO TO 100
ZMIN=PMTS(JSIDE)/(W(1)*W(2))
ZMAX=1.-PMTS(3-JSIDE)/(W(1)*W(2))
IF(ZMIN.GE.ZMAX) GO TO 100
C..WARNING: VERY IMPORTANT THE ORDER OF IFL AND IFLN IN ZFRAGS
c.. fraction of momentum acc. to the fragmentation fct.
Z=ZFRAGS(IFL(JSIDE),IFLN,PTS,ZMIN,ZMAX,leading)
PPTCL(3,I)=0.5*(Z*W(JSIDE)-PMTS(JSIDE)/
*(Z*W(JSIDE)))*(-1.)**(JSIDE+1)
PPTCL(4,I)=0.5*(Z*W(JSIDE)+PMTS(JSIDE)/(Z*W(JSIDE)))
IDCAY(I)=0
IF(.NOT.(JSIDE.EQ.1)) GO TO 282
IDENTR(NPR)=IDENT(I)
PPTR(1,NPR)=PPTCL(1,I)
PPTR(2,NPR)=PPTCL(2,I)
PPTR(3,NPR)=PPTCL(3,I)
PPTR(4,NPR)=PPTCL(4,I)
PPTR(5,NPR)=PPTCL(5,I)
282 IF(.NOT.(JSIDE.EQ.2)) GO TO 283
IDENTL(NPL)=IDENT(I)
PPTL(1,NPL)=PPTCL(1,I)
PPTL(2,NPL)=PPTCL(2,I)
PPTL(3,NPL)=PPTCL(3,I)
PPTL(4,NPL)=PPTCL(4,I)
PPTL(5,NPL)=PPTCL(5,I)
283 IF(DIQBR) GO TO 210
IFL(JSIDE)=-IFLN
PX1(JSIDE)=-PX2
PY1(JSIDE)=-PY2
GO TO 220
C NEW DIQUARK CREATION
210 ID1=IABS(IFLL)
ID2=IABS(IFLN)
IF(ID1.LE.ID2) GO TO 215
ISWAP=ID1
ID1=ID2
ID1=ISWAP
215 IFL(JSIDE)=ISIGN(1000*ID1+100*ID2,IFLL)
PX1L(JSIDE)=PX1L(JSIDE)+PXL-PX2
PY1L(JSIDE)=PY1L(JSIDE)+PYL-PY2
PX1(JSIDE)=PX1L(JSIDE)
PY1(JSIDE)=PY1L(JSIDE)
220 W(1)=W(1)-PPTCL(4,I)-PPTCL(3,I)
W(2)=W(2)-PPTCL(4,I)+PPTCL(3,I)
SPINT=.TRUE.
KSPIN=2
IF(MOD(IFL(1),100).EQ.0.AND.MOD(IFL(2),100).EQ.0) GO TO 240
IDB=IDPARS(IFL(1),IFL(2),SPINT,KSPIN)
AMB=AMASS(IDB)+dmass
GO TO 211
240 IFCN=1
IF(ranf(0).GT.0.5) IFCN=2
IFLC1=-IFCN
IF(IFL(1).GT.0) IFLC1=IFCN
IFLC2=-IFLC1 IKH1=IDPARS(IFL(1),IFLC1,SPINT,KSPIN)
IKH2=IDPARS(IFL(2),IFLC2,SPINT,KSPIN)
AMB=AMASS(IKH1)+AMASS(IKH2)+DMASs
211 P1X=PX1(1)+PX1(2)
P1Y=PY1(1)+PY1(2)
PT12=P1X**2+P1Y**2
W12=W(1)*W(2)
AMS2=W12-PT12
IF(AMS2.LT.AMB**2) GO TO 100
SPINT=.TRUE.
KSPIN=1
IF(MOD(IFL(1),100).EQ.0.AND.MOD(IFL(2),100).EQ.0) GO TO 231
IDR=IDPARS(IFL(1),IFL(2),SPINT,KSPIN)
WEND=(AMASS(IDR)+DMAS)**2
GO TO 232
231 IKHR1=IDPARS(IFL(1),IFLC1,SPINT,KSPIN)
IKHR2=IDPARS(IFL(2),IFLC2,SPINT,KSPIN)
WEND=(AMASS(IKHR1)+AMASS(IKHR2)+DMAS)**2
232 SPINT=.FALSE.
KSPIN=0
IF(W(1)*W(2).GE.WEND) GO TO 130
GO TO 230
225 P1X=PX1(1)+PX1(2)
P1Y=PY1(1)+PY1(2)
PT12=P1X**2+P1Y**2
W12=W(1)*W(2)
AMS2=W12-PT12
C LAST BREAK OF STRING
230 NPTCL=I
AMC=SQRT(AMS2)
EC=(W(1)+W(2))/2.0
VC(1)=P1X/EC
VC(2)=P1Y/EC
VC(3)=(W(1)-W(2))/(2.0*EC)
NIN1=NPTCL+1
c.. the last break of the string will be done in clustr
CALL CLUSTR(IFL(1),IFL(2),AMC,ilead)
IF(LRET) GO TO 100
NFIN1=NPTCL
CALL LORTR(VC,NIN1,NFIN1,BACK)
NPR=NPR+1
NPL=NPL+1
IF(NPR.GT.MXPTCL.OR.NPL.GT.MXPTCL) GO TO 9999
c..the hadron from the left and the right side of the string
c..are copied to the final pptcl array
IDENTL(NPL)=IDENT(NFIN1)
PPTL(1,NPL)=PPTCL(1,NFIN1)
PPTL(2,NPL)=PPTCL(2,NFIN1)
PPTL(3,NPL)=PPTCL(3,NFIN1)
PPTL(4,NPL)=PPTCL(4,NFIN1)
PPTL(5,NPL)=PPTCL(5,NFIN1)
IDENTR(NPR)=IDENT(NIN1)
PPTR(1,NPR)=PPTCL(1,NIN1)
PPTR(2,NPR)=PPTCL(2,NIN1)
PPTR(3,NPR)=PPTCL(3,NIN1)
PPTR(4,NPR)=PPTCL(4,NIN1)
PPTR(5,NPR)=PPTCL(5,NIN1)
JJ=NFIX
DO 284 J=1,NPR
JJ=JJ+1
IDENT(JJ)=IDENTR(J)
PPTCL(1,JJ)=PPTR(1,J)
PPTCL(2,JJ)=PPTR(2,J)
PPTCL(3,JJ)=PPTR(3,J)
PPTCL(4,JJ)=PPTR(4,J)
PPTCL(5,JJ)=PPTR(5,J)
284 CONTINUE
JJ=NFIX+NPR
DO 285 J=1,NPL
JJ=JJ+1 K=NPL-J+1
IDENT(JJ)=IDENTL(K)
PPTCL(1,JJ)=PPTL(1,K)
PPTCL(2,JJ)=PPTL(2,K)
PPTCL(3,JJ)=PPTL(3,K)
PPTCL(4,JJ)=PPTL(4,K)
PPTCL(5,JJ)=PPTL(5,K)
285 CONTINUE
N1=NFIX+1
N2=NFIX+NPR+NPL-1
c.. we choose the LUND scheme
consti=.false.
IF(CONSTI) THEN
C------------------------------------------------------C
C----- CONSTITUENT TIME ----------------C
C------------------------------------------------------C
DO 1286 J=N1,N2
P3S=0.
ES=0.
DO 1287 L=N1,J
P3S=P3S+PPTCL(3,L)
1287 ES=ES+PPTCL(4,L)
c.. TI is the formation time of the particle
c.. ZI is the z coordinate
TI=(AMSTR-2.*P3S)/(2.*XAP)
ZI=(AMSTR-2.*ES)/(2.*XAP)
IF(J.NE.N2) GO TO 1288
TII=TI
ZII=ZI
1288 PPTCL(6,J)=0.
PPTCL(7,J)=0.
PPTCL(8,J)=ZI
PPTCL(9,J)=TI
1286 CONTINUE
C
PPTCL(6,N2+1)=0.
PPTCL(7,N2+1)=0.
PPTCL(8,N2+1)=ZII
PPTCL(9,N2+1)=TII
C
GO TO 1253
ENDIF
C------------------------------------------------------C
C----- INSIDE-OUTSIDE TIME (LUND) ----------------C
C------------------------------------------------------C
DO 286 J=N1,NPTCL
P3S=0.
ES=0.
NJ=J-1
IF(NJ.EQ.0) GO TO 289
DO 287 L=N1,NJ
P3S=P3S+PPTCL(3,L)
287 ES=ES+PPTCL(4,L)
c.. TI is the formation time of the particle
c.. ZI is the z coordinate
289 TI=(AMSTR-2.*P3S+PPTCL(4,J)-PPTCL(3,J))/(2.*XAP)
ZI=(AMSTR-2.*ES-PPTCL(4,J)+PPTCL(3,J))/(2.*XAP)
PPTCL(6,J)=0.
PPTCL(7,J)=0.
PPTCL(8,J)=ZI
PPTCL(9,J)=TI
286 CONTINUE
1253 RETURN
c.. warning if to many hadrons are produced in string
c.. increase the particle arrays to avoid this
9999 WRITE(6,9998) I
9998 FORMAT(//10X,40H...STOP IN STRING..NPTCL TOO HIGH NPTCL=,I5)
STOP END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
SUBROUTINE CLUSTR(IFL1,IFL2,AMCTR,ilead)
c
cinput amctr : stringmass,
cinput ifl1 : leading quark (or diquark) along $+Z$ axis
cinput ifl2 : 2nd leading quark (or diquark)
cinput ilead : $2-ilead=$ number of leading (di-)quarks
c
c output : produced particles via common block ({\tt pptcl})
c
C HADRONS PRODUCTION BY MEANS CLUSTER BREAKING
C WITH QUARK AND ANTIQUARK OR QUARK AND DIQUARK OR DIQUARK AND
C ANTIDIQUARK IFL1 AND IFL2 ON ENDS.
c Only the final 2 particles are created in {\tt clustr}!
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
implicit real*8 (a-h,o-z)
implicit integer (i-n)
COMMON/COMTRY/ NTRIES
PARAMETER(MXPTCL=200)
COMMON/PARTCL/ PPTCL(9,MXPTCL),nptcl,IDENT(MXPTCL),IDCAY(MXPTCL)
include 'comstr.f'
COMMON/COLRET/ LRET
LOGICAL LRET
DIMENSION IFL(2),U(3)
LOGICAL SPINT
real*8 valint(1)
common /values/ valint
include 'options.f'
ctp060202 to avoid warnings with gfortran compilation
logical ctp060202
ctp060202=.false.
if(ctp060202)write(*,*)ilead
ctp060202 end
c.. strange and charm suppression in clustr (see string)
prbs=ctparam(6)
prbc=ctparam(7)
c..the diquark-suppresion parameter is reduced for small
c string masses (finite size effect) see A. Jahns diploma thesis
if(amctr.le.2.8d0)then
pbarc=0.d0
elseif (amctr.le.5d0)then
pbarc=((amctr-2.8d0)/2.2d0)**3*ctparam(8)
else
pbarc=ctparam(8)
endif
dmas=ctparam(9)
dmass=ctparam(10)
C
C PRBS STRANGENESS SUPPRESSION PARAMETER
C PRBC CHARM SUPPRESSION PARAMETER
PU=1./(2.+PRBS+PRBC)
C
NFIX=NPTCL
NREP=0
LRET=.FALSE.
100 I=NFIX
IF(NREP.LT.NTRIES) GO TO 101
LRET=.TRUE. RETURN
101 CONTINUE
KSPIN=0
IFL(1)=IFL1
IFL(2)=IFL2
SPINT=.FALSE.
I=I+2
IF(I.GT.MXPTCL) GO TO 9999
C CHOOSE SIDE OF BREAK
JSIDE=1
C IF ANY IFL IS A DIQUARK
IF(MOD(IFL(1),100).EQ.0.OR.MOD(IFL(2),100).EQ.0) GO TO 150
C IFL(1) AND IFL(2) ARE QUARKS
C SELECT Q,QBARPAIR OR QQ,QQBAR PAIR
DRND=ranf(0)
IF(DRND.LT.PBARC.and.valint(1).eq.0.d0) GO TO 140
C Q,QBAR PAIR
IFLN=ISIGN(IFLAV(PU,PRBS),-IFL(JSIDE))
GO TO 200
C QQ,QQBAR PAIR
140 IQ1=IFLAV(PU,PRBS)
IQ2=IFLAV(PU,PRBS)
IF(IQ1.LE.IQ2) GO TO 145
ISWAP=IQ1
IQ1=IQ2
IQ2=ISWAP
145 IFQQ=1000*IQ1+100*IQ2
IFLN=ISIGN(IFQQ,IFL(JSIDE))
GO TO 200
C IFL(1) OR IFL(2) IS DIQUARK
C Q,QBAR PAIR
150 IPSIGN=IFL(JSIDE)
IF(MOD(IFL(JSIDE),100).EQ.0) GO TO 130
IPSIGN=-IFL(JSIDE)
130 IFLN=ISIGN(IFLAV(PU,PRBS),IPSIGN)
C IDENTS AND MASSES OF PARTICLES
200 continue
c..quark-flip included (to describe some phi data)
if(CTParam(5).gt.ranf(0).and.mod(ifln,100).ne.0.and.
& mod(ifl(jside),100).ne.0.and.mod(ifl(3-jside),100).ne.0)then
c quark-flip
IDENT(I-1)=IDPARC(IFL(JSIDE),IFL(3-JSIDE),SPINT,KSPIN)
IDENT(I)=IDPARC(-IFLN,IFLN,SPINT,KSPIN)
else
IDENT(I-1)=IDPARC(IFL(JSIDE),IFLN,SPINT,KSPIN)
IDENT(I)=IDPARC(IFL(3-JSIDE),-IFLN,SPINT,KSPIN)
end if
c..for special bbar-b annihilation reactions for conservation
c of total quantum numbers
if(valint(1).ne.0.d0)then
ifq1=int(valint(1)/10.)
ifq2=-mod(int(valint(1)),10)
if(isign(1,ifln).eq.isign(1,ifq1))then
IDENT(I-1)=IDPARC(IFL(JSIDE),ifq1,SPINT,KSPIN)
IDENT(I)=IDPARC(IFL(3-JSIDE),ifq2,SPINT,KSPIN)
else
IDENT(I-1)=IDPARC(IFL(JSIDE),ifq2,SPINT,KSPIN)
IDENT(I)=IDPARC(IFL(3-JSIDE),ifq1,SPINT,KSPIN)
endif
endif
PPTCL(5,I-1)=AMASS(IDENT(I-1))
PPTCL(5,I)=AMASS(IDENT(I))
C IF TOO LOW MASS,START ALL OVER (i.e. goto 100)
DEMAS=0.15
IF(IFLN.LT.3) DEMAS=0.
c IF(AMCTR.GT.PPTCL(5,I-1)+PPTCL(5,I)+DEMAS) GO TO 102
IF(AMCTR.GT.PPTCL(5,I-1)+PPTCL(5,I)+DEMAS) then
if(mod(ifl1,100).eq.0.or.mod(ifl2,100).eq.0)goto 102
c..maximum kinetic energy cutoff for meson-clustr:c..we want a lot of energy in the particle mass in this last break
IF(AMCTR-PPTCL(5,I-1)-PPTCL(5,I).lt.ctparam(43)) GO TO 102
endif
NREP=NREP+1
c.. 100 starts all over
GO TO 100
102 CONTINUE
c.. isotropic px py pz distribution
PA=DBLPCM(AMCTR,PPTCL(5,I-1),PPTCL(5,I))
U(3)=1.-2.*ranf(0)
PHI=2.*PI*ranf(0)
ST=SQRT(1.-U(3)**2)
U(1)=ST*COS(PHI)
U(2)=ST*SIN(PHI)
PPTCL(1,I-1)=PA*U(1)
PPTCL(1,I)=-(PA*U(1))
PPTCL(2,I-1)=PA*U(2)
PPTCL(2,I)=-(PA*U(2))
PPTCL(3,I-1)=PA*U(3)
PPTCL(3,I)=-(PA*U(3))
PA2=PA**2
PPTCL(4,I-1)=SQRT(PA2+PPTCL(5,I-1)**2)
PPTCL(4,I)=SQRT(PA2+PPTCL(5,I)**2)
IDCAY(I-1)=0
IDCAY(I)=0
NPTCL=I
c..forward/backward distribution in clustr for baryons
c..(no pt in the last string break!)
c..pt for the baryon comes from parton kick in the excitation
if(abs(ident(i)).ge.1000.or.abs(ident(i-1)).ge.1000)then
PPTCL(1,I-1)=0.d0
PPTCL(1,I)=0.d0
PPTCL(2,I-1)=0.d0
PPTCL(2,I)=0.d0
PPTCL(3,I-1)=PA
PPTCL(3,I)=-PA
PA2=PA**2
PPTCL(4,I-1)=SQRT(PA2+PPTCL(5,I-1)**2)
PPTCL(4,I)=SQRT(PA2+PPTCL(5,I)**2)
IDCAY(I-1)=0
IDCAY(I)=0
NPTCL=I
endif
c..if baryon number=+-1, force the (anti-)baryon in positive
c..z-direction (just pick the right hemisphere)
if(ctoption(29).gt.0)then
if( (iabs(ident(i)/1000).ne.0.and.iabs(ident(i-1)/1000).eq.0
& .and.pptcl(3,i).lt.0.d0).or.
& (iabs(ident(i-1)/1000).ne.0.and.iabs(ident(i)/1000).eq.0
& .and.pptcl(3,i-1).lt.0.d0)) then
pptcl(3,i) =-pptcl(3,i)
pptcl(3,i-1)=-pptcl(3,i-1)
endif
endif
RETURN
c.. particle array to small warning:
9999 WRITE(6,9998) I
9998 FORMAT(//10X,40H...STOP IN CLUSTR..NPTCL TOO HIGH NPTCL=,I5)
STOP
END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
integer FUNCTION IFLAV(PU,PRBS)
c
cinput PU : 1-{\tt PU}= up (down, resp.) probabilitycinput PRBS : Strange quark suppression
c
c output : {\tt iflav}: flavor of created quark
c
c Returns quark flavor acc. to suppression prob's:
c 1=up, 2=down, 3=strange, 4=charm
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
implicit real*8 (a-h,o-z)
implicit integer (i-n)
C
RNDOM=ranf(0)
C
IF(RNDOM.GT.PU) GO TO 1
c..create up quark
IFLAV=1
RETURN
1 IF(RNDOM.GT.2.0*PU) GO TO 2
c..create down quark
IFLAV=2
RETURN
2 IF(RNDOM.GT.PU*(2.0+PRBS)) GO TO 3
c..create strange quark
IFLAV=3
RETURN
c..create charm quark
3 IFLAV=4
RETURN
END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
real*8 FUNCTION ZFRAGS(IFL,IFLN,PT2,ZMIN,ZMAX,leading)
c
cinput IFL : ID of existing quark
cinput IFLN : ID of newly created quark
cinput PT2 : $p_t$ of newly created hadron
cinput ZMIN : lowest allowed longitudinal momentum fraction
cinput ZMAX : highest allowed longitudinal momentum fraction
cinput leading : flag for leading particle
c
coutput ZFRAGS : longitudinal momentum fraction of created hadron
c
c According to the fragmentation function(s), longitudinal momentum
c is assigned to the hadron.
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
implicit real*8 (a-h,o-z)
implicit integer (i-n)
LOGICAL leading
include 'options.f'
COMMON/INPRNT/ ITDKY,ITLIS
PARAMETER(ALFT=0.5,ARHO=0.5,APHI=0.,APSI=-2.)
PARAMETER(AN=-0.5,ALA=-0.75,ALAC=-1.75)
PARAMETER(AKSI=-1.0,AUSC=-2.0,AUCC=-2.0)
c.. cto 21 chooses the fragmentation fct.
if ((.not.leading).or.(ifln.eq.3)) then
affm=ctparam(47)
bffm=ctparam(48)
5108 ZFRAGS=ZMIN+ranf(0)*(ZMAX-ZMIN)
YF=((1d0-ZFRAGS)**bffm*(bffm+1)*affm+1-affm)/3d0
ctp060926 if(yf.gt.1d0.or.yf.lt.0d0)then
ctp060926 write(6,*)'ZFRAGS: wrong norm:',yf,zmin,zmax,zfrags
ctp060926 end if
IF(ranf(0).LE.YF) RETURN
GO TO 5108
endif
c.. GAUSSIAN fragmentation fct.
if(CTOption(21).eq.0)then
affm=CTParam(36)
bffm=CTParam(37)
c.. suppress low momentum particles
deltaz=zmax-zmin
zmin=zmin+deltaz*0.25d0
108 ZFRAGS=ZMIN+ranf(0)*(ZMAX-ZMIN)
c.. gaussian distr.
yf=exp(-(((zfrags-bffm)**2)/(2.*affm**2)))
c.. field-feynmann fragmentation fct.
c YF=((1d0-ZFRAGS)**bffm*(bffm+1)*affm+1-affm)/3d0
ctp060926 if(yf.gt.1d0.or.yf.lt.0d0)then
ctp060926 write(6,*)'ZFRAGS: wrong norm:',yf,zmin,zmax,zfrags
ctp060926 end if
c return
IF(ranf(0).LE.YF) RETURN
GO TO 108
else if(CTOption(21).eq.1)then
c..lund-fragmentation fct.
1008 ZFRAGS=ZMIN+ranf(0)*(ZMAX-ZMIN)
YF=(1-zfrags)*exp(-(0.7*pt2/zfrags))/3d0
ctp060926 if(yf.gt.1d0.or.yf.lt.0d0)then
ctp060926 write(6,*)'ZFRAGS: wrong norm:',yf,zmin,zmax,zfrags
ctp060926 end if
IF(ranf(0).LE.YF) RETURN
GO TO 1008
else if(CTOption(21).eq.2)then
c.. kaidalov's fragmentation fct.
ID1=IABS(IFL)
ID2=IABS(IFLN)
IF(MOD(ID2,100).EQ.0) GO TO 15
GO TO(1,2,3,4),ID2
C UU-TRAJECTORY
1 ZFRAGS=ZMIN+ranf(0)*(ZMAX-ZMIN)
YF=(1.0-ZFRAGS)**(ALFT-ARHO)
IF(ranf(0).LE.YF) RETURN
GO TO 1
C DD-TRAJECTORY
2 ZFRAGS=ZMIN+ranf(0)*(ZMAX-ZMIN)
YF=(1.-ZFRAGS)**(ALFT-ARHO)
IF(ranf(0).LE.YF) RETURN
GO TO 2
C SS-TRAJECTORY
3 ZFRAGS=ZMIN+ranf(0)*(ZMAX-ZMIN)
YF=(1.-ZFRAGS)**(ALFT-APHI)
IF(ranf(0).LE.YF) RETURN
GO TO 3
C CC-TRAJECTORY
4 ZFRAGS=ZMIN+ranf(0)*(ZMAX-ZMIN)
YF=(1.-ZFRAGS)**(ALFT-APSI)
IF(ranf(0).LE.YF) RETURN
GO TO 4
C
15 Continue
CALL FLAVOR(ID2,IFL2,IFL3,IFL1,ISPIN)
IF((IFL2.EQ.1.AND.IFL3.EQ.1)) GO TO 16
IF((IFL2.EQ.1.AND.IFL3.EQ.2)) GO TO 17
IF((IFL2.EQ.1.AND.IFL3.EQ.3)) GO TO 18
IF((IFL2.EQ.1.AND.IFL3.EQ.4)) GO TO 19
IF((IFL2.EQ.2.AND.IFL3.EQ.2)) GO TO 20
IF((IFL2.EQ.2.AND.IFL3.EQ.3)) GO TO 21
IF((IFL2.EQ.2.AND.IFL3.EQ.4)) GO TO 22
IF((IFL2.EQ.3.AND.IFL3.EQ.3)) GO TO 23
IF((IFL2.EQ.3.AND.IFL3.EQ.4)) GO TO 24
IF((IFL2.EQ.4.AND.IFL3.EQ.4)) GO TO 25
C UUUU-TRAJECTORY16 ZFRAGS=ZMIN+ranf(0)*(ZMAX-ZMIN)
YF=(1.-ZFRAGS)**(ALFT-(2.*AN-ARHO))
IF(ranf(0).LE.YF) RETURN
GO TO 16
C UDUD-TRAJECTORY
17 ZFRAGS=ZMIN+ranf(0)*(ZMAX-ZMIN)
YF=(1.-ZFRAGS)**(ALFT-(2.*AN-ARHO))
IF(ranf(0).LE.YF) RETURN
GO TO 17
C USUS-TRAJECTORY
18 ZFRAGS=ZMIN+ranf(0)*(ZMAX-ZMIN)
YF=(1.-ZFRAGS)**(ALFT-(2.*ALA-ARHO))
IF(ranf(0).LE.YF) RETURN
GO TO 18
C UCUC-TRAJECTORY
19 ZFRAGS=ZMIN+ranf(0)*(ZMAX-ZMIN)
YF=(1.-ZFRAGS)**(ALFT-(2.*ALAC-ARHO))
IF(ranf(0).LE.YF) RETURN
GO TO 19
C DDDD-TRAJECTORY
20 ZFRAGS=ZMIN+ranf(0)*(ZMAX-ZMIN)
YF=(1.-ZFRAGS)**(ALFT-(2.*AN-ARHO))
IF(ranf(0).LE.YF) RETURN
GO TO 16
C DSDS-TRAJECTORY
21 ZFRAGS=ZMIN+ranf(0)*(ZMAX-ZMIN)
YF=(1.-ZFRAGS)**(ALFT-(2.*ALA-ARHO))
IF(ranf(0).LE.YF) RETURN
GO TO 21
C DCDC-TRAJECTORY
22 ZFRAGS=ZMIN+ranf(0)*(ZMAX-ZMIN)
YF=(1.-ZFRAGS)**(ALFT-(2.*ALAC-ARHO))
IF(ranf(0).LE.YF) RETURN
GO TO 22
C SSSS-TRAJECTORY
23 ZFRAGS=ZMIN+ranf(0)*(ZMAX-ZMIN)
YF=(1.-ZFRAGS)**(ALFT-(2.*AKSI-ARHO))
IF(ranf(0).LE.YF) RETURN
GO TO 23
C SCSC-TRAJECTORY
24 ZFRAGS=ZMIN+ranf(0)*(ZMAX-ZMIN)
YF=(1.-ZFRAGS)**(ALFT-(2.*AUSC-ARHO))
IF(ranf(0).LE.YF) RETURN
GO TO 24
C CCCC-BARYON
25 ZFRAGS=ZMIN+ranf(0)*(ZMAX-ZMIN)
YF=(1.-ZFRAGS)**(ALFT-(2.*AUCC-ARHO))
IF(ranf(0).LE.YF) RETURN
GO TO 25
else
write(6,*)'string.f: cto(21)=',ctoption(21),' not valid'
stop
end if
END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
SUBROUTINE LORTR(V,NIN,NFIN,BACK)
c
cinput V : boost velocity (3-vector)
cinput NIN : lower boundary in the pptcl array
cinput NFIN : upper boundary in the pptcl array
cinput back : inversion flag for transformation
c
c output : via common block ({\tt pptcl})
c
c Performs a Lorentz-boost of a part of the pptcl array
c (from {\tt nin} to {\tt nfin})
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
implicit real*8 (a-h,o-z) implicit integer (i-n)
PARAMETER(MXPTCL=200)
COMMON/PARTCL/ PPTCL(9,MXPTCL),nptcl,IDENT(MXPTCL),IDCAY(MXPTCL)
DIMENSION V(3),VD(3)
LOGICAL BACK
L=1
IF(BACK) L=-1
DO 3 I=1,3
3 VD(I)=V(I)
VVD=VD(1)*VD(1)+VD(2)*VD(2)+VD(3)*VD(3)
GAD=1.D0/DSQRT(DABS(1.D0-VVD))
GA=GAD
DO 100 J=NIN,NFIN
VP=V(1)*PPTCL(1,J)+V(2)*PPTCL(2,J)+V(3)*PPTCL(3,J)
GAVP=GA*(GA*VP/(1.+GA)-FLOAT(L)*PPTCL(4,J))
PPTCL(1,J)=PPTCL(1,J)+GAVP*V(1)
PPTCL(2,J)=PPTCL(2,J)+GAVP*V(2)
PPTCL(3,J)=PPTCL(3,J)+GAVP*V(3)
PMAS=PPTCL(5,J)
PPTCL(4,J)=SQRT(PPTCL(1,J)**2+PPTCL(2,J)**2+PPTCL(3,J)**2+
+PMAS**2)
100 CONTINUE
RETURN
END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
BLOCK DATA D2
c
c Initial values for several common blocks
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
implicit real*8 (a-h,o-z)
implicit integer (i-n)
C
C INPUT QUARK DISTRIBUTION PARAMETERS
C
include 'comstr.f'
COMMON/KAPPA/ XAP
COMMON/CONSTI/ CONSTI
LOGICAL CONSTI
DATA CONSTI/.FALSE./
C
C
C DATA FOR COSPAR
c this parameter controls the min. energy for ending string fragmentation
c default DATA DMAS/0.35/
c now: ctparam(9) and ctparam(10)
c DATA DMAS/1./
c data dmass/0.25/
C
C
C INPUT STRING TENSION (FERMI/GEV)
DATA XAP/1.0/
C
C INPUT FRAGMENTATION PARAMETERS:
C
DATA PARQLS/0.5/,PARRS/0.9/
C PRBS STRANGENESS SUPPRESSION PARAMETER
C PRBC CHARM SUPPRESSION PARAMETER
C PBARC AND PBARS BARYON PAIRS SUPPRESION PARAMETERS
C DATA PRBS/0.3/,PRBC/0.0005/,PBARC/0.09/,PBARS/0.09/
c these are now ctparam(6),ctparam(7) and ctparam(8)
c DATA PRBS/0.33/,PRBC/0.0/,PBARC/0.09/,PBARS/0.09/
C
C THESE PARAMETERS ARE IMPORTANT FOR RESONANCE PRODUCTION
DATA PJSPNC/.50/
DATA PJSPNS/.50/
c DATA PJSPNC/1./
c DATA PJSPNS/1./C DATA PMIX1C/.25,.25,.5,.5,.5,1./,PMIX2C/.5,.5,1.,0.,0.,1./
c DATA PMIX1S/.25,.25,.5,.5,.5,1./,PMIX2S/.5,.5,1.,0.,0.,1./
c QUARK MIXING PARAMETERS
cspl-0795 these parameters assign the pure u/ubar,d/dbar,s/sbar states
c (e.g. 110,220,330) to the physical particles according to the su(3)
c quark model. The flavor mixing angles are chosen according to quadratic
c Gell-Mann-Okubo mass formula (Review of Particle Properties,
c Phys. Rev D50 (1994) 1319). For the scalar mesons this formula is not
c applicable. We assume an ideal mixing angle (tan(theta)=1/sqrt(2)).
c
c pseudoscalar: theta=-10 deg
c vector : theta= 39 deg
c pseudovector: theta= 51 deg
c tensor : theta= 28 deg
c
c PMIX1C/S(IF1,JSPIN) gives the quark content of the heavy isosinglet:
c e.g. the eta' (330) is 25% u/ubar 25% d/dbar and 50% s/sbar
c PMIX2C/S(IF1,JSPIN) chooses for nonstrange qqbars between the isosinglet
c the and the isotriplet
c quark content of the heavy singlett
c DATA PMIX1C/.25,.25,.5, ! eta'
c & 0.,0.,1., ! phi
c & 0.,0.,1., ! f_0
c & .04,.04,.92, ! f_1'
c & .01,.01,.98/ ! f_2'
c DATA PMIX2C/.5 ,.5 ,1.,.5,.5,1.,.5,.5,1.,.5 ,.5 ,1. ,.5 ,.5 , 1./
c DATA PMIX1S/.25,.25,.5,0.,0.,1.,0.,0.,1.,.04,.04,.92,.01,.01,.98/,
c & PMIX2S/.5 ,.5 ,1.,.5,.5,1.,.5,.5,1.,.5 ,.5 ,1. ,.5 ,.5 , 1./
C
C SIGQTS IS IMPORTANT PARAMETER FOR PRODUCED HADRON TRANSVERSE MOMENTA
C
c DATA SIGQTS/0.65/ -> ctparam(42)
END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
BLOCK DATA D4
c
cinput () : NONE
c
coutput () : via common blocks
c
c Initial values for several common blocks
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
implicit real*8 (a-h,o-z)
implicit integer (i-n)
COMMON/COMTRY/ NTRIES
C
C
C DATA COMTRY
DATA NTRIES/1000/
C
END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
integer FUNCTION IDPARS(IFL01,IFL02,SPINT,IR)
c
cinput IFL01 : ID of (anti-)quarks/di-quarks
cinput IFL02 : ID of (anti-)quarks/di-quarks
cinput SPINT : flag for spin assignment
cinput IR : Determines particle spin
c
c output : Quark code of the hadron
c
C CONSTRUCT MESON FROM QUARK AND ANTIQUARK WITH FLAVORS IFL01,IFL02
C OR CONSTRUCT BARYON FROM DIQUARK AND QUARK OR ANTIDIQUARK AND
C ANTIQUARK WITH FLAVORS IFL01,IFL02.
c THE MESON MULTIPLETT IS CHOSEN ACC. TO SUPPRESSION PARAM'S:
c parm gives the probability for different meson multiplets according
c to spin degeneracy and average mass ratios
c spin-parity 0- : 1- : 0+ : 1+ : 2+ = parm(1):parm(2)...:parm(5)c If SPINT=.t., IR will be used to assign particle spin
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
implicit real*8 (a-h,o-z)
implicit integer (i-n)
LOGICAL SPINT
include 'options.f'
include 'comstr.f'
C
IFL1=IFL01
IFL2=IFL02
C CONSTRUCT MESON WITH ACCOUNT FLAVOR MIXING
IF(MOD(IFL1,100).EQ.0) GO TO 420
IF(MOD(IFL2,100).EQ.0) GO TO 425
c
call getbran(parm,1,njspin,psdum,1,njspin,jspin)
jspin=jspin-1
IF(SPINT.AND.IR.EQ.2) JSPIN=0
IF(SPINT.AND.IR.EQ.1) JSPIN=1
ID1=IFL1
ID2=IFL2
IF(ID1+ID2.NE.0) GO TO 400
C
ID1=IABS(ID1)
IF(ID1.GE.4) GO TO 401
RND=ranf(0)
ID1=INT(PMIX1S(ID1,JSPIN+1)+RND)+INT(PMIX2S(ID1,JSPIN+1)+RND)+1
401 ID2=-ID1
C
400 IF(IABS(ID1).LE.IABS(ID2)) GO TO 410
ISAVE=ID1
ID1=ID2
ID2=ISAVE
410 IDHAD=ISIGN(100*IABS(ID1)+10*IABS(ID2)+JSPIN,ID1)
GO TO 470
C CONSTRUCT BARYON IDENT
420 ID3=ISIGN(MOD(IFL1/100,10),IFL1)
ID2=IFL1/1000
ID1=IFL2
GO TO 430
425 ID3=ISIGN(MOD(IFL2/100,10),IFL2)
ID2=IFL2/1000
ID1=IFL1
430 IF(IABS(ID1).LE.IABS(ID2)) GO TO 431
ISWAP=ID1
ID1=ID2
ID2=ISWAP
431 IF(IABS(ID2).LE.IABS(ID3)) GO TO 432
ISWAP=ID2
ID2=ID3
ID3=ISWAP
432 IF(IABS(ID1).LE.IABS(ID2)) GO TO 440
ISWAP=ID1
ID1=ID2
ID2=ISWAP
440 JSPIN=1
IF(ID1.EQ.ID2.AND.ID2.EQ.ID3) GO TO 450
JSPIN=INT(ranf(0)+PJSPNS)
IF(SPINT.AND.IR.EQ.2) JSPIN=0
IF(SPINT.AND.IR.EQ.1) JSPIN=1
450 IF(JSPIN.EQ.1.OR.ID1.EQ.ID2.OR.ID2.EQ.ID3) GO TO 460
DRND=ranf(0)
IF(DRND.GT.PJSPNS) GO TO 460
ISWAP=ID1 ID1=ID2
ID2=ISWAP
460 IDHAD=1000*IABS(ID1)+100*IABS(ID2)+10*IABS(ID3)+JSPIN
IDHAD=ISIGN(IDHAD,IFL1)
470 IDPARS=IDHAD
RETURN
END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
integer FUNCTION IDPARC(IFL01,IFL02,SPINT,IR)
c
cinput IFL01 : ID of (anti-)quarks/di-quarks
cinput IFL02 : ID of (anti-)quarks/di-quarks
cinput SPINT : flag for spin assignment
cinput IR : Determines particle spin
c
c output : Quark code of the hadron
c
C CONSTRUCT MESON FROM QUARK AND ANTIQUARK WITH FLAVORS IFL01,IFL02
C OR CONSTRUCT BARYON FROM DIQUARK AND QUARK OR ANTIDIQUARK AND
C ANTIQUARK WITH FLAVORS IFL01,IFL02.
c THE MESON MULTIPLETT IS CHOSEN ACC. TO SUPPRESSION PARAM'S:
c parm gives the probability for different meson multiplets according
c to spin degeneracy and average mass ratios
c spin-parity 0- : 1- : 0+ : 1+ : 2+ = parm(1):parm(2)...:parm(5)
c If SPINT=.t., IR will be used to assign particle spin
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
implicit real*8 (a-h,o-z)
implicit integer (i-n)
LOGICAL SPINT
include 'options.f'
include 'comstr.f'
C
IFL1=IFL01
IFL2=IFL02
C CONSTRUCT MESON WITH ACCOUNT FLAVOR MIXING
IF(MOD(IFL1,100).EQ.0) GO TO 420
IF(MOD(IFL2,100).EQ.0) GO TO 425
c.. choose multiplett by its probability
call getbran(parm,1,njspin,dummy,1,njspin,jspin)
jspin=jspin-1
IF(SPINT.AND.IR.EQ.2) JSPIN=0
IF(SPINT.AND.IR.EQ.1) JSPIN=1
ID1=IFL1
ID2=IFL2
IF(ID1+ID2.NE.0) GO TO 400
C
ID1=IABS(ID1)
IF(ID1.GE.4) GO TO 401
c.. singlet mixing acc. to mixing angles
RND=ranf(0)
ID1=INT(PMIX1C(ID1,JSPIN+1)+RND)+INT(PMIX2C(ID1,JSPIN+1)+RND)+1
401 ID2=-ID1
C
400 IF(IABS(ID1).LE.IABS(ID2)) GO TO 410
ISAVE=ID1
ID1=ID2
ID2=ISAVE
410 IDHAD=ISIGN(100*IABS(ID1)+10*IABS(ID2)+JSPIN,ID1)
GO TO 470
C CONSTRUCT BARYON IDENT
420 ID3=ISIGN(MOD(IFL1/100,10),IFL1)
ID2=IFL1/1000
ID1=IFL2
GO TO 430
425 ID3=ISIGN(MOD(IFL2/100,10),IFL2) ID2=IFL2/1000
ID1=IFL1
430 IF(IABS(ID1).LE.IABS(ID2)) GO TO 431
ISWAP=ID1
ID1=ID2
ID2=ISWAP
431 IF(IABS(ID2).LE.IABS(ID3)) GO TO 432
ISWAP=ID2
ID2=ID3
ID3=ISWAP
432 IF(IABS(ID1).LE.IABS(ID2)) GO TO 440
ISWAP=ID1
ID1=ID2
ID2=ISWAP
440 JSPIN=1
IF(ID1.EQ.ID2.AND.ID2.EQ.ID3) GO TO 450
JSPIN=INT(ranf(0)+PJSPNC)
IF(SPINT.AND.IR.EQ.2) JSPIN=0
IF(SPINT.AND.IR.EQ.1) JSPIN=1
450 IF(JSPIN.EQ.1.OR.ID1.EQ.ID2.OR.ID2.EQ.ID3) GO TO 460
DRND=ranf(0)
IF(DRND.GT.PJSPNC) GO TO 460
ISWAP=ID1
ID1=ID2
ID2=ISWAP
460 IDHAD=1000*IABS(ID1)+100*IABS(ID2)+10*IABS(ID3)+JSPIN
IDHAD=ISIGN(IDHAD,IFL1)
470 IDPARC=IDHAD
RETURN
END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
real*8 FUNCTION AMASS(ID)
c
c
cinput ID : Quark code
c
c
C THIS FUNCTION RETURNS THE MASS OF THE PARTICLE WITH
C IDENT CODE ID. (QUARK-BASED IDENT CODE)
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
implicit real*8 (a-h,o-z)
implicit integer (i-n)
include 'comres.f'
include 'comstr.f'
real*8 mmin,mmax,m0,mminit, massit, widit
integer isoit
dimension amq(4)
include 'options.f'
c.. quark masses (u,d,s,c)
DATA AMq/.15,.15,.45,1.6/
c fraction of baryonresonances
bresfrac=ctparam(11)
CALL FLAVOR(ID,IFL1,IFL2,IFL3,JSPIN)
idabs=iabs(id)
c get quark masses
if (idabs.le.4) then
amass=amq(idabs)
return
endif
c get diquark masses
IF(ID.NE.0.AND.MOD(ID,100).EQ.0) then
AMASS=AMq(IABS(IFL1))+AMq(IABS(IFL2))
return
endif
c get hadron masses
c get baryon masses
c (anti-)nucleon ?
if ((idabs.eq.1120).or.(idabs.eq.1220)) then
if (ranf(0).lt.bresfrac)then
c.. N*
amass=getmass(mresmax,1)
return
else
c.. N
amass=0.938
return
endif
endif
c (anti-)delta ?
if ((idabs.eq.1111).or.(idabs.eq.1121)
& .or.(idabs.eq.1221).or.(idabs.eq.2221)) then
if (ranf(0).lt.bresfrac)then
c.. D*
amass=getmass(mresmax,2)
return
else
c.. Delta(1232)
m0=massit(mindel)
w0=widit(mindel)
c get meson mass accord. to breit wigner distr.
mmin=mminit(mindel)
mmax=m0+3d0*w0
cdh call getmas(m0,w0,mindel,isoit(mindel),mmin,mmax,-1,amass)
call getmas(m0,w0,mindel,isoit(mindel),mmin,mmax,-1.d0,amass)
return
endif
endif
c (other baryons)
if(idabs-1000.ge.0) then
call id2ityp(id,0.d0,itypin,iz2)
c..check range and avoid double counting for explicitely treated resonances
if ((abs(itypin).ge.minlam.and.abs(itypin).le.maxcas).and.
& (idabs.ne.1231.and.idabs.ne.1131.and.idabs.ne.2231.and.
& idabs.ne.1331.and.idabs.ne.2331).and.
& (ctoption(31).eq.1)) then
c.. generate also non-groundstate ityps for strange baryons
call probitypres(itypin,ityp)
else
ityp=itypin
endif
amass=massit(iabs(ityp))
return
endif
c mesons
if(idabs.le.330+njspin-1) then
call id2ityp(id,0.d0,ityp,iz2)
m0=massit(iabs(ityp))
w0=widit(iabs(ityp))
c get meson mass accord. to breit wigner distr.
mmin=max(mminit(iabs(ityp)),m0-3d0*w0)
mmax=m0+3d0*w0
cdh call getmas(m0,w0,ityp,iz2,mmin,mmax,-1,amass)
call getmas(m0,w0,ityp,iz2,mmin,mmax,-1.d0,amass)
return
endif
write(6,*)'! amass: no mass for part.id:',id
return
END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
real*8 FUNCTION DBLPCM(A,B,C)c
cinput A : Mass of particle A
cinput B : Mass of particle B
cinput C : Mass of particle C
c
c
c In the rest frame of a particle of mass A, decaying into particles
c of masses B and C, {\tt dblpcm} returns the momenta of the outgoing
c particles.
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
implicit real*8 (a-h,o-z)
implicit integer (i-n)
DA=A
DB=B
DC=C
DVAL=(DA**2-DB**2-DC**2)**2-(2.D0*DB*DC)**2
DBLPCM=0.
IF(DVAL.GT.0.D0)DBLPCM=DSQRT(DVAL)/(2.D0*DA)
return
END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine ityp2id(ityp,iz2,ifa,ifb)
c
cinput ityp : UrQMD particle ID
cinput iz2 : UrQMD $2\cdot I_3$
c
coutput ifa : quarkcode (diquark)
coutput ifb : quarkcode (quark)
c
c returns quark id from uqmd-ityp and isospin z-component (times 2)
c
c quark ids are: 1 up, 2 down, 3 strange, 4 charm
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
implicit none
include 'comres.f'
integer itypabs,ityp,iz2,ifa,ifb,sf
integer t3,if(3)
itypabs=iabs(ityp)
t3=iz2*isign(1,ityp)
c nucleons ?
if (itypabs.lt.minmes) then
sf=strres(itypabs)
c uuu
if(t3.eq.3) then
if(1)=1
if(2)=1
if(3)=1
call mquarks(if,ifa,ifb)
ifa=ifa*isign(1,ityp)
ifb=ifb*isign(1,ityp)
return
endif
c uus
if(t3.eq.2) then
if(1)=1
if(2)=1
if(3)=3
call mquarks(if,ifa,ifb)
ifa=ifa*isign(1,ityp)
ifb=ifb*isign(1,ityp)
return
endif
c uud
if((t3.eq.1).and.(sf.eq.0)) then
if(1)=1
if(2)=1 if(3)=2
call mquarks(if,ifa,ifb)
ifa=ifa*isign(1,ityp)
ifb=ifb*isign(1,ityp)
return
endif
c uss
if((t3.eq.1).and.(sf.eq.2)) then
if(1)=1
if(2)=3
if(3)=3
call mquarks(if,ifa,ifb)
ifa=ifa*isign(1,ityp)
ifb=ifb*isign(1,ityp)
return
endif
c uds
if((t3.eq.0).and.(sf.eq.1)) then
if(1)=1
if(2)=2
if(3)=3
call mquarks(if,ifa,ifb)
ifa=ifa*isign(1,ityp)
ifb=ifb*isign(1,ityp)
return
endif
c sss
if((t3.eq.0).and.(sf.eq.3)) then
if(1)=3
if(2)=3
if(3)=3
call mquarks(if,ifa,ifb)
ifa=ifa*isign(1,ityp)
ifb=ifb*isign(1,ityp)
return
endif
c udd
if((t3.eq.-1).and.(sf.eq.0)) then
if(1)=1
if(2)=2
if(3)=2
call mquarks(if,ifa,ifb)
ifa=ifa*isign(1,ityp)
ifb=ifb*isign(1,ityp)
return
endif
c dss
if((t3.eq.-1).and.(sf.eq.2)) then
if(1)=2
if(2)=3
if(3)=3
call mquarks(if,ifa,ifb)
ifa=ifa*isign(1,ityp)
ifb=ifb*isign(1,ityp)
return
endif
c dds
if (t3.eq.-2) then
if(1)=2
if(2)=2
if(3)=3
call mquarks(if,ifa,ifb)
ifa=ifa*isign(1,ityp)
ifb=ifb*isign(1,ityp)
return
endif
c ddd
if (t3.eq.-3) then
if(1)=2
if(2)=2 if(3)=2
call mquarks(if,ifa,ifb)
ifa=ifa*isign(1,ityp)
ifb=ifb*isign(1,ityp)
return
endif
endif
c-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
c bosons
if (itypabs.ge.minmes) then
sf=strmes(itypabs)
c d ubar
if(t3.eq.-2) then
ifa=2
ifb=-1
return
endif
c d sbar
if((t3.eq.-1).and.(sf*isign(1,ityp).eq.-1)) then
ifa=2
ifb=-3
return
endif
c s ubar
if((t3.eq.-1).and.(sf*isign(1,ityp).eq.+1)) then
ifa=3
ifb=-2
return
endif
c q qbar
if(t3.eq.0) then
c u ubar
c all neutral triplet states:
if ((itypabs.eq.minmes+1) .or. !pi0
& (itypabs.eq.minmes+4) .or.
& (itypabs.eq.minmes+11) .or.
& (itypabs.eq.minmes+14) .or.
& (itypabs.eq.minmes+18) .or.
& (itypabs.eq.minmes+22) .or.
& (itypabs.eq.minmes+26) .or.
& (itypabs.eq.minmes+30) .or.
c the gamma gets also a quark content
& (itypabs.eq.minmes)) then
ifa=1
ifb=-1
return
endif
c d dbar
c light singlet states (generally less strange quark content)
if ((itypabs.eq.minmes+2) .or. !eta
& (itypabs.eq.minmes+3) .or.
& (itypabs.eq.minmes+5) .or.
& (itypabs.eq.minmes+15) .or.
& (itypabs.eq.minmes+19) .or.
& (itypabs.eq.minmes+23) .or.
& (itypabs.eq.minmes+27) .or.
& (itypabs.eq.minmes+31)) then
ifa=2
ifb=-2
return
endif
c s sbar
if ((itypabs.eq.minmes+7) .or. !eta' (958)
& (itypabs.eq.minmes+9) .or. !phi (1020)
& (itypabs.eq.minmes+12) .or. !f_0 (980)
& (itypabs.eq.minmes+16) .or. !f_1 (1510)
& (itypabs.eq.minmes+20) .or. !f_2'(1525)
& (itypabs.eq.minmes+24) .or. !f_2'(1525)
& (itypabs.eq.minmes+28) .or. !f_2'(1525) & (itypabs.eq.minmes+32))then
ifa=3
ifb=-3
return
endif
endif
c u sbar
if((t3.eq.1).and.(sf*isign(1,ityp).eq.-1)) then
ifa=1
ifb=-3
return
endif
c s dbar
if((t3.eq.1).and.(sf*isign(1,ityp).eq.+1)) then
ifa=3
ifb=-1
return
endif
c u dbar
if(t3.eq.2) then
ifa=1
ifb=-2
return
endif
endif
c in any other case we will do a mesonic string
c.. and a warning !
ctp060926 write(*,*)'! ityp2id: ityp',ityp,',iz2',iz2,
ctp060926 & ' can not be converted into id. Please check.'
ifa=1
ifb=-1
RETURN
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine mquarks(if,ifa,ifb)
c
cinput if : single quarks (array)
c
coutput ifa : diquark
coutput ifb : quark
c
c this routine adds randomly the single quarks of the
c baryon to become a diquark and a quark.
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
implicit none
include 'options.f'
real*8 ranf
integer if(0:2),ifa,ifb,ir
integer i,prod,inons
ir=int(3.*ranf(0))
ifa=1000*if(ir)+100*if(mod(ir+1,3))
ifb=if(mod(ir+2,3))
c.. check if heavy quark clusters are switched off (cto 37=0)
if (CTOption(37).eq.0) return
c.. not switched off -> clusters strange quarks to
c.. diquark molecule, i.e. keep ss-diquark together
c. are there 2 strange quarks
prod=if(0)*if(1)*if(2)
if (prod.eq.9.or.prod.eq.18.or.prod.eq.27) then
c. find the non-strange quark inons
inons=0 do i=0,2
if(if(i).lt.3) inons=i
enddo
ifb=if(inons)
ifa=1000*if(mod(inons+1,3))+100*if(mod(inons+2,3))
endif
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine id2ityp(id,mass,it,iz)
c
cinput id : quarkcode
cinput mass : mass of resonance
c
coutput: it : UrQMD particle ID
coutput: iz : $2\cdot I_3$ of particle
c
c returns UrQMD ID ({\tt it}) and isospin z-component (times 2) ({\tt IZ})
c from quark id
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
implicit none
include 'comres.f'
include 'options.f'
integer IT,IZ,id,idabs,i,j,k,jspin,idloc,whichres
real*8 mass,pardelt,mminit,dm,dmold,massit
integer iit,irun
c.. extract quark id's and spin
IDABS=IABS(ID)
I=IDABS/1000
J=MOD(IDABS/100,10)
K=MOD(IDABS/10,10)
JSPIN=MOD(IDABS,10)
c single diquarks are not popular in urqmd:
IF(ID.NE.0.AND.MOD(ID,100).EQ.0) GO TO 222
c single quarks are not popular in urqmd:
IF(J.EQ.0) GO TO 222
c mesons:
if(I.EQ.0) then
c calculate isospin from quark content:
IZ=0
if(j.lt.3) IZ=3-2*j
if(k.lt.3) IZ=IZ+(2*k-3)
IZ=isign(1,id)*IZ
idloc=idabs-jspin
if((idloc.eq.110).or.(idloc.eq.120))then
c..triplett (e.g. pion)
IT= mlt2it(jspin*4+1) ! minmes+1
return
else if(idloc.eq.220)then
c..1st singlett (e.g. eta)
IT= mlt2it(jspin*4+3) !minmes+2
return
else if(idloc.eq.130.or.idloc.eq.230)then
c..strange doublett (e.g. (anti-)kaon)
IT=isign(mlt2it(jspin*4+2),id)
return
else if(idloc.eq.330)then
c..2nd singlett (e.g. eta')
IT=mlt2it(jspin*4+4) !minmes+7
return
else
goto 222
endif
else
c.. baryons:
c calculate isospin from quark content:
IZ=0
if(i.lt.3) IZ=3-2*i
if(j.lt.3) IZ=IZ+3-2*j
if(k.lt.3) IZ=IZ+3-2*k
IZ=isign(1,id)*IZ
c spin 1/2 baryons (all nucleon-resonances are treated here!)
if(jspin.eq.0)then
c (anti-)nucleons and resonances
if(idabs.eq.1120.or.idabs.eq.1220)then
c mass below parnuc -> nucleon (ground state), mass above parnuc ->N*
if(mass.lt.mminit(minnuc+1))then
IT=isign(minnuc,id)
return
else
IT=isign(whichres(dble(mass),1),id)
return
endif
else if(idabs.eq.2130)then
c lambda
iit=minlam
if (mass.gt.1d0)then
dmold=1d30
do irun = minlam,maxlam
dm=abs(massit(irun)-mass)
if (dm.le.dmold)then
dmold=dm
iit=irun
end if
end do
end if
IT=isign(iit,id)
return
else if(idabs.eq.1230.or.idabs.eq.1130.or.idabs.eq.2230)then
c sigma
iit=minsig
if (mass.gt.1d0)then
dmold=1d30
do irun = minsig,maxsig
dm=abs(massit(irun)-mass)
if (dm.le.dmold)then
dmold=dm
iit=irun
end if
end do
end if
IT=isign(iit,id)
return
else if(idabs.eq.1330.or.idabs.eq.2330)then
c cascade
iit=mincas
if (mass.gt.1d0)then
dmold=1d30
do irun = mincas,maxcas
dm=abs(massit(irun)-mass)
if (dm.le.dmold)then
dmold=dm
iit=irun
end if
end do
endif
IT=isign(iit,id)
return
else
goto 222
endif
c spin 3/2 baryons (all delta-resonances are treated here!) else if(jspin.eq.1)then
if(idabs.eq.1111.or.idabs.eq.1121.or.
& idabs.eq.1221.or.idabs.eq.2221)then
c if mass below pardelt -> Delta1232, otherwise: Delta-resonance
pardelt=1.45
if(mass.lt.pardelt)then
c delta 1232
IT=isign(whichres(dble(mass),0),id)
return
else
c delta resonance
IT=isign(whichres(dble(mass),2),id)
return
endif
else if(idabs.eq.1231.or.idabs.eq.1131.or.idabs.eq.2231)then
c sigma*
IT=isign(minsig+1,id)
return
else if(idabs.eq.1331.or.idabs.eq.2331)then
c cascade*
IT=isign(mincas+1,id)
return
else if(idabs.eq.3331)then
c omega
IT=isign(minome,id)
return
else
goto 222
endif
c higher spin baryons include here:
else
goto 222
endif
endif
222 continue
ctp060926 write(6,*)'! ID=',id,' can not be converted into ityp'
ctp060926 write(6,*)'I=',i,'J=',j,'K=',k,'spin=',jspin
RETURN
END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine probitypres(itypin,itypout)
c
cinput itypin : ityp of groundstate
c
coutput: itypout : ityp of groundstate plus resonances
c
c returns a new ityp which also includes resonce ityps, this is
c necessary to include hyperon resonces as long as no proper getmass
c for hyperons existst.
c the probabilities are chooses according to a exp(delta m/b) distribution.
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
implicit none
include 'comres.f'
integer itypin,itypout,ss,i,itypinn
real*8 prob(1:200),probres,norm,deltam,massit
ss=itypin/abs(itypin)
itypinn=abs(itypin)
c.. the lambda
if (itypinn.ge.minlam.and.itypinn.le.maxlam)then
c get probabilities and norm
norm=0d0
do i=minlam+1,maxlam
deltam=abs(massit(i)-massit(minlam)) prob(i)=probres(deltam,massit(minlam),i)
norm=norm+prob(i)
end do
do i=minlam+1,maxlam
prob(i)=prob(i)/norm
end do
call findityp(prob,itypout,minlam+1,maxlam)
itypout=itypout*ss
return
end if
c.. the sigma
if (itypinn.ge.minsig.and.itypinn.le.maxsig)then
c get probabilities and norm
norm=0d0
do i=minsig+1,maxsig
deltam=abs(massit(i)-massit(minsig))
prob(i)=probres(deltam,massit(minsig),i)
c... do not generate masses for minsig+1 they are treated ecplicitely
if (i.eq.minsig+1) prob(i)=0d0
norm=norm+prob(i)
end do
do i=minsig+1,maxsig
prob(i)=prob(i)/norm
end do
call findityp(prob,itypout,minsig+1,maxsig)
itypout=itypout*ss
return
end if
c.. the cascades
if (itypinn.ge.mincas.and.itypinn.le.maxcas)then
c get probabilities and norm
norm=0d0
do i=mincas+1,maxcas
deltam=abs(massit(i)-massit(mincas))
prob(i)=probres(deltam,massit(mincas),i)
c... do not generate masses for mincas+1 they are treated ecplicitely
if (i.eq.mincas+1) prob(i)=0d0
norm=norm+prob(i)
end do
do i=mincas+1,maxcas
prob(i)=prob(i)/norm
end do
call findityp(prob,itypout,mincas+1,maxcas)
itypout=itypout*ss
return
end if
write(*,*)'itypin, itypout',itypin,itypout
stop 'Error in probitypres!'
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
real*8 function probres(dm,minmass,it)
c
cinput dm : mass difference to groundstate
cinput minmass : mass of groundstate
cinput it : ityp
c
coutput: probres: unnormalized probability for this state
c
c returns probabilty for a higher mass state according to
c exponential distribution with degeneracy
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
implicit none
real*8 dm,minmass,j,g,T
integer it, getspin
c.. assume temperature T=170 MeV (from statistical model, e.g. becattini,heinz)
T=0.170d0
c.. get spin*2
j=1d0*getspin(it,1)
if (j.lt.0) then
c.. take spin into account via J= m^2 and deg. g=2j+1 (regge theory)
j=2d0*(minmass+dm)**2
end if
g=(j+1d0)
probres=g*exp(-(dm/T))
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine findityp (p,it,mini,maxi)
c
cinput p : array with normalized probabilities
cinput mini : lowest index
cinput maxi : largest index
c
coutput: it: ityp according to probability
c
c returns a new ityp according to probabilties defined in p
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
implicit none
real*8 p(1:200),y,ranf
integer it,mini,maxi,ix
it=mini
ix=0
y=0d0
1 continue
ix=int(mini+int((maxi-mini+1)*ranf(0)))
y=ranf(0)
if (y.lt.p(ix)) then
it=ix
return
end if
goto 1
end
C---------------------------------------------------------------------------
C THE END