File indexing completed on 2024-04-06 12:14:05
0001
0002 subroutine utresc(iret)
0003
0004
0005
0006
0007 include 'epos.inc'
0008 include 'epos.incems'
0009 double precision p1,esoll,ppp,seedp,psoll,pt1soll,pt2soll
0010 dimension p1(5),p2(4),p0(5,mamx+mamx),pini(mxptl)
0011 logical force,nolead(mxptl),lspec(mxptl),lim
0012 data scalmean/0./scalevt/0./
0013 save scalmean,scalevt
0014 call utpri('utresc',ish,ishini,4)
0015
0016 errlim=0.005 !max(0.001,1./engy)
0017 if(iLHC.eq.1)errlim=max(0.00005,0.5/engy)
0018
0019 iret=0
0020 nptlpt=iabs(maproj)+iabs(matarg)
0021 call ranfgt(seedp) !not to change the seed ...
0022 if(nptl.le.nptlpt) goto 9999
0023
0024 if(ish.ge.8)then
0025 call alistf('list before boost&')
0026 endif
0027 esoll=0.d0
0028 psoll=0.d0
0029 p1(1)=0.d0
0030 p1(2)=0.d0
0031 p1(3)=0.d0
0032 p1(4)=0.d0
0033 p2(3)=0.d0
0034 p2(4)=0.d0
0035 ipmax=4
0036 imin=nptlpt+1
0037 if(iappl.eq.1)then
0038 imin=1
0039 ipmax=2
0040 if(iLHC.eq.1)ipmax=0
0041 endif
0042
0043 do i=1,nptlpt
0044 nolead(i)=.false.
0045 do j=1,5
0046 p0(j,i)=dble(pptl(j,i))
0047 enddo
0048
0049
0050 do j=ipmax+1,4
0051 p1(j)=p1(j)+dble(pptl(j,i))
0052 enddo
0053
0054
0055 if(mod(istptl(i),10).eq.0)then
0056 do j=ipmax+1,4
0057 p2(j)=p2(j)+dble(pptl(j,i))
0058 enddo
0059 endif
0060 enddo
0061
0062 do i=nptlpt+1,nptl
0063 if(mod(istptl(i),10).eq.0)then
0064
0065 if(iLHC.eq.1.and.pptl(4,i).gt.engy*0.51)then
0066 if(ish.ge.1)write(ifch,*)'Energy of particle too high !'
0067 & ,i,ityptl(i),pptl(4,i)
0068
0069 if(ityptl(i).eq.48.or.ityptl(i).eq.58 !diffractive resonance
0070 & .or.ityptl(i).eq.47.or.ityptl(i).eq.57)then !active spectators
0071 pptl(4,i)=0.5*engy
0072 amt=sqrt(pptl(1,i)**2+pptl(2,i)**2+pptl(5,i)**2)
0073 pptl(3,i)=(pptl(4,i)+amt)*(pptl(4,i)-amt)
0074 if(pptl(3,i).gt.0.)then
0075 pptl(3,i)=sqrt(pptl(3,i))
0076 else
0077 iret=1
0078 endif
0079 else
0080 iret=1
0081 endif
0082 if(iret.eq.1)then
0083 if(ish.ge.1)write(ifch,*)'Warning in utresc: redo event...'
0084
0085 goto 9999
0086 endif
0087 endif
0088
0089 do j=1,ipmax
0090 p1(j)=p1(j)+dble(pptl(j,i))
0091 enddo
0092
0093 do j=ipmax+1,4
0094 p2(j)=p2(j)+dble(pptl(j,i))
0095 enddo
0096 lspec(i)=.false.
0097 if(((ityptl(i).eq.45.or.ityptl(i).eq.48).and.maproj.ge.100)
0098 & .or.((ityptl(i).eq.55.or.ityptl(i).eq.58).and.matarg.ge.100))
0099 & lspec(i)=.true.
0100 if((abs(pptl(3,i)/pnullx).le.0.9
0101 & .and.abs(pptl(3,i)).gt.pptl(5,i)).or.lspec(i))then
0102 nolead(i)=.true.
0103
0104 else
0105 nolead(i)=.false.
0106
0107 endif
0108 endif
0109 enddo
0110 psoll=max(dble(pnullx),abs(p1(3)))
0111
0112 if(iappl.eq.1)then
0113 diff=abs(p2(3)-p1(3))
0114 scal=p2(4)/p1(4)
0115 if(abs(scal-1.).le.errlim.and.abs(diff/psoll).lt.errlim
0116 & .and.(iLHC.eq.0.or.
0117 & (abs(p2(1)).lt.errlim.and.abs(p2(2)).lt.errlim)))then
0118 if(ish.ge.4)
0119 & write (ifch,'(2x,a,2g14.6)') 'Energy OK: ',scal,abs(diff/psoll)
0120 goto 9999
0121 else
0122 diff=0.
0123 scal=1.
0124 endif
0125 endif
0126
0127 ppp=(p1(4)+p1(3))*(p1(4)-p1(3))-p1(2)*p1(2)-p1(1)*p1(1)
0128 if(ppp.gt.0.d0)then
0129 p1(5)=sqrt(ppp)
0130 else
0131 iret=1
0132 write(ifch,*)'p1',p1(1),p1(2),p1(3),p1(4),ppp
0133 if(ish.ge.2)write (ifch,*) 'Problem in utresc (1): redo event'
0134 write (ifmt,*) 'Problem in utresc (1): redo event'
0135
0136 goto 9999
0137 endif
0138 esoll=p1(5)
0139 if(ish.ge.4) write (ifch,'(a,5g14.6)') 'boost-vector: ',p1
0140
0141
0142
0143 pmax=0.d0
0144 npart=0
0145 pt1soll=0.d0
0146 pt2soll=0.d0
0147 do i=imin,nptl
0148 if(mod(istptl(i),10).le.1)then
0149 call utlob4(1,p1(1),p1(2),p1(3),p1(4),p1(5)
0150 $ ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i))
0151 if(mod(istptl(i),10).eq.0.and.i.gt.nptlpt)then
0152 npart=npart+1
0153 pt1soll=pt1soll+dble(pptl(1,i))
0154 pt2soll=pt2soll+dble(pptl(2,i))
0155 endif
0156 endif
0157 if(i.gt.nptlpt.and.nolead(i))pmax=max(pmax,abs(pptl(3,i)))
0158 enddo
0159
0160 if(ish.ge.6)then
0161 call alistf('list after boost&')
0162 endif
0163
0164 if(ish.ge.5)write(ifch,'(a)')'--------rescale momenta----------'
0165
0166 if(iLHC.eq.1)then
0167
0168 if(ish.ge.6)write(ifch,*)'ptsoll:',pt1soll,pt2soll,npart
0169 pt1soll=pt1soll/dble(npart)
0170 pt2soll=pt2soll/dble(npart)
0171 do i=nptlpt+1,nptl
0172 if(mod(istptl(i),10).eq.0)then
0173 pptl(1,i)=pptl(1,i)-sngl(pt1soll)
0174 pptl(2,i)=pptl(2,i)-sngl(pt2soll)
0175 pptl(4,i)=sqrt(pptl(1,i)**2+pptl(2,i)**2
0176 & +pptl(3,i)**2+pptl(5,i)**2)
0177 endif
0178 enddo
0179
0180 endif
0181
0182 if(ish.ge.6)write(ifch,*)'esoll,psoll,pmax:',esoll,psoll,pmax
0183
0184
0185
0186 scal=1.
0187 diff0=0.
0188
0189 ferr=0.05
0190 force=.false.
0191 npart=nptl-imin+1
0192 do ipass=1,300
0193 sum=0.
0194 sum3=0.
0195 difft=diff0
0196 ndif=0
0197 nfirst=int(rangen()*float(npart)) !start list randomly
0198 do i=0,npart-1
0199 j=imin+i+nfirst
0200 if(j.gt.nptl)j=imin+i+nfirst-npart
0201 if(mod(istptl(j),10).eq.0)then
0202
0203 if(nolead(j))then
0204
0205
0206 if(scal.eq.1..and.abs(diff0).lt.1.e-6)then
0207 ndif=ndif+1
0208 pini(j)=pptl(3,j)
0209 else
0210 pptl3new=0.
0211 if( force .or.(
0212 & ityptl(j)/10.eq.4.or.ityptl(j)/10.eq.5
0213 & ))then !try just remnant first
0214 ndif=ndif+1
0215 diff=sign(min(0.3*abs(pini(j)),
0216 & rangen()*abs(difft)),difft)
0217 pptl3new=scal*(pptl(3,j)-diff)
0218
0219
0220 if(abs(pptl3new).lt.pmax)then
0221
0222 if(abs(pptl3new-pini(j)).lt.ferr*abs(pini(j))
0223 & .or.(lspec(j).and.abs(pptl3new).lt.abs(0.8*pini(j))))then
0224
0225 difft=difft-diff
0226 pptl(3,j)=scal*(pptl(3,j)-diff)
0227 pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2
0228 * +pptl(3,j)**2+pptl(5,j)**2)
0229 endif
0230 endif
0231 endif
0232 endif
0233 endif
0234
0235 sum=sum+pptl(4,j)
0236 sum3=sum3+pptl(3,j)
0237 endif
0238 enddo
0239
0240 diff=sum3
0241 scal=sngl(esoll)/(sum-diff)
0242 if(ish.ge.6)write(ifch,*)
0243 $ 'ipass,scal,diff/psoll,e,pz,ndif,f:'
0244 $ ,ipass,scal,diff/psoll,sum,sum3,ndif,force
0245 if(abs(scal-1.).le.errlim.and.abs(diff/psoll).lt.errlim)
0246 $ goto 300
0247 if(ndif.gt.0.and.(force.or.ipass.lt.150))then
0248
0249 diff0=diff
0250 elseif(abs(scal-1.).le.1e-2.and.abs(diff/psoll).lt.5e-2
0251 & .and.iLHC.eq.0)then
0252 goto 300
0253 elseif(force)then
0254 if(ish.ge.2)
0255 $ write(ifmt,*)'Warning in utresc: no more allowed particle'
0256 goto 302
0257 else
0258 force=.true.
0259 ferr=0.1
0260 diff=0.
0261 endif
0262 enddo
0263 302 if(iLHC.eq.1)then
0264 lim=.not.(abs(scal-1.).le.errlim.and.abs(diff/psoll).lt.errlim)
0265 else
0266 lim=abs(scal)+abs(diff/psoll).gt.2.5
0267 endif
0268
0269 if(ish.ge.1)then
0270 call utmsg('utrescl')
0271 write(ifch,*)'***** scal=',scal,diff/psoll,lim
0272 call utmsgf
0273 endif
0274
0275 if(lim)then
0276 if(ish.ge.1)then
0277 write(ifmt,*)'Warning in utresc !'
0278 write(ifch,'(a,i10,d25.15)')'redo EPOS event ...'
0279 & ,nint(seedj),seedc
0280 endif
0281 iret=1
0282 goto 9999
0283 endif
0284
0285
0286
0287 300 continue
0288
0289 do i=1,nptl
0290 if(i.le.nptlpt)then
0291 do j=1,5
0292 pptl(j,i)=p0(j,i)
0293 enddo
0294 else
0295 if(mod(istptl(i),10).le.1)then
0296 call utlob4(-1,p1(1),p1(2),p1(3),p1(4),p1(5)
0297 $ ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i))
0298 endif
0299 endif
0300 enddo
0301
0302 if(ish.ge.4)call alist('list after rescaling&',1,nptl)
0303
0304 9999 continue
0305 if(ish.ge.2)then
0306 scalevt=scalevt+1.
0307 scalmean=scalmean+scal
0308 write(ifch,*)' average rescaling factor: ',scalmean
0309 & /scalevt
0310 endif
0311 call ranfst(seedp) ! ... after this subroutine
0312 call utprix('utresc',ish,ishini,4)
0313
0314 end
0315
0316
0317 subroutine utghost(iret)
0318
0319
0320
0321 include 'epos.inc'
0322 include 'epos.incems'
0323 double precision seedp
0324 call utpri('ughost',ish,ishini,4)
0325
0326 iret=0
0327 nptlpt=iabs(maproj)+iabs(matarg)
0328 if(iappl.eq.6.or.iappl.eq.8)nptlpt=3 ! ee or DIS
0329 call ranfgt(seedp) !not to change the seed ...
0330 if(nptl.le.nptlpt) goto 9999
0331
0332 if(ish.ge.5)write(ifch,'(a)')'---------mark ghosts---------'
0333
0334
0335
0336 do j=nptlpt+1,nptl
0337 if(istptl(j).le.1.and.pptl(4,j).gt.0.d0)then
0338 if(iLHC.eq.0.or.mod(abs(idptl(j)),10).le.1)then !for LHC tune don't fix mass of resonnances (to keep width)
0339 amass=pptl(5,j)
0340 call idmass(idptl(j),amass)
0341 if(abs(idptl(j)).gt.100.and.
0342 & abs(pptl(5,j)-amass).gt.0.01*amass)then
0343 if(ish.ge.5)write(ifch,*)'wrong particle mass',j,idptl(j)
0344 & ,pptl(5,j),amass
0345 amass=pptl(5,j)
0346 call idres(idptl(j),amass,idr,iadj)
0347 if(idr.ne.0)then
0348 pptl(5,j)=amass
0349 idptl(j)=idr
0350 else
0351 call idmass(idptl(j),amass)
0352 pptl(5,j)=amass
0353 endif
0354 call idtau(idptl(j),pptl(4,j),pptl(5,j),taugm)
0355 tivptl(2,j)=tivptl(1,j)+taugm*(-alog(rangen()))
0356 else
0357 pptl(5,j)=amass
0358 endif
0359 endif
0360 if(abs((pptl(4,j)+pptl(3,j))*(pptl(4,j)-pptl(3,j))
0361 $ -pptl(2,j)**2-pptl(1,j)**2-pptl(5,j)**2).gt.0.3
0362 $ .and.abs(1.-abs(pptl(3,j))/pptl(4,j)).gt.0.01)then
0363 !print*,'ghost',ityptl(j),idptl(j)
0364 if(ish.ge.1)write(ifmt,*)'ghost:',j,idptl(j),ityptl(j)
0365 if(ish.ge.5)then
0366 write(ifch,'(a,$)')'ghost:'
0367 call alistc("&",j,j)
0368 endif
0369 ityptl(j)=100+ityptl(j)/10
0370 elseif(irescl.ge.1)then
0371 c ensure that all particles are really on-shell
0372 pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2
0373 * +pptl(3,j)**2+pptl(5,j)**2)
0374 endif
0375 elseif(mod(istptl(j),10).eq.0)then
0376 c if not droplet with fusion
0377 if(istptl(j).ne.10.or.iorsdf.ne.3)then
0378 if(ish.ge.1)then
0379 write(ifmt,*)'Lost particle (E=0)'
0380 write(ifch,*)'Lost particle (E=0) :'
0381 call alistc("utghost&",j,j)
0382 endif
0383 istptl(j)=istptl(j)+2
0384 endif
0385 endif
0386 enddo
0387
0388 if(ish.ge.5)write(ifch,'(a)')'---------treat ghosts---------'
0389
0390 c treat ghosts
0391 c ------------
0392 ifirst=1
0393 scal=1.
0394 pfif=0.
0395 efif=0.
0396 ntry=0
0397 132 nfif=0
0398 psum=0
0399 esum=0.
0400 ntry=ntry+1
0401 do j=nptlpt+1,nptl
0402 if(mod(istptl(j),10).eq.0)then
0403 if(ityptl(j).gt.100)then
0404 nfif=nfif+1
0405 if(ifirst.eq.1)then
0406 pfif=pfif+pptl(3,j)
0407 if(pptl(4,j).gt.0.)efif=efif+pptl(4,j)
0408 endif
0409 if(irescl.ge.1) then
0410 if(ifirst.gt.1)then
0411 if(pptl(4,j).gt.0.)then
0412 Einv=1./pptl(4,j)
0413 amt=1.-(pptl(5,j)*Einv)**2+(pptl(1,j)*Einv)**2
0414 $ +(pptl(2,j)*Einv)**2
0415 else
0416 amt=-1.
0417 endif
0418 if(amt.gt.0.)then
0419 pptl(3,j)=sign(pptl(4,j),pptl(3,j))*sqrt(amt)
0420 else
0421 y=(rangen()+rangen()+rangen()+rangen()-2.)/2.*yhaha
0422 y=sign(abs(y),pptl(3,j))
0423 pptl(3,j)
0424 $ =sqrt(pptl(5,j)**2+pptl(1,j)**2
0425 $ +pptl(2,j)**2)*sinh(y)
0426 pptl(4,j)
0427 $ =sqrt(pptl(5,j)**2+pptl(1,j)**2
0428 $ +pptl(2,j)**2)*cosh(y)
0429 efif=efif+pptl(4,j)
0430 endif
0431 ifirst=0
0432 else
0433 c do k=1,3
0434 do k=3,3
0435 pptl(k,j)=pptl(k,j)*scal
0436 enddo
0437 pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2
0438 * +pptl(5,j)**2)
0439 endif
0440 endif
0441 psum=psum+pptl(3,j)
0442 esum=esum+pptl(4,j)
0443 if(ish.ge.5)
0444 $ write (ifch,*) 'nrevt,psum,esum,pfif,efif,nfif,scal'
0445 $ ,nrevt,psum,esum,pfif,efif,nfif,scal
0446 endif
0447 endif
0448 enddo
0449 if ( ish.ge.5 ) write (ifch,*) 'tot',nfif,efif,pfif,esum,psum
0450
0451
0452 if(nfif.gt.5.or.(esum.gt.0.05*engy.and.nfif.ne.1))then
0453 if(ifirst.eq.0)then
0454 do j=nptlpt+1,nptl
0455 if ( ityptl(j).ge.101 .and. ityptl(j).le.105 )then
0456 if((psum-pfif)*(1.-scal).ge.0)
0457 & pptl(3,j)=pptl(3,j)-(psum-pfif)/nfif
0458 endif
0459 enddo
0460 else
0461 ifirst=2
0462 goto 132
0463 endif
0464 scal=efif/esum
0465 if ( ish.ge.5 ) write (ifch,*) 'scal',scal
0466 if ( abs(scal-1.) .gt. 0.05 ) then
0467 if(ntry.le.1000)then
0468 goto 132
0469 else
0470 iret=1
0471 if(ish.ge.2)write (ifch,*) 'Problem in utghost : redo event'
0472 if(ish.ge.1)write (ifmt,*) 'Problem in utghost : redo event'
0473 goto 9999
0474 endif
0475 endif
0476 else
0477 do j=nptlpt+1,nptl
0478 if ( ityptl(j).ge.101 .and. ityptl(j).le.105 )then
0479 pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2
0480 * +pptl(5,j)**2)
0481 endif
0482 enddo
0483 endif
0484
0485 if(ish.ge.5)write(ifch,'(a)')'---------Check Ghost list---------'
0486
0487 c Check Ghost list
0488
0489 if(ish.ge.5)then
0490 do j=nptlpt+1,nptl
0491 if(mod(istptl(j),10).eq.0)then
0492 if(ityptl(j).le.105.and.ityptl(j).ge.101)then
0493 write(ifch,'(a,$)')'ghost:'
0494 call alistc("&",j,j)
0495 endif
0496 endif
0497 enddo
0498 endif
0499
0500
0501 9999 continue
0502 call ranfst(seedp) ! ... after this subroutine
0503 call utprix('ughost',ish,ishini,4)
0504
0505 end
0506
0507 c-----------------------------------------------------------------------
0508 subroutine utrsph(iret)
0509 c-----------------------------------------------------------------------
0510 c if irescl=1 and ispherio=1 rescaling is done for particle used by
0511 c spherio as initial condition.
0512 c-----------------------------------------------------------------------
0513 include 'epos.inc'
0514 include 'epos.incems'
0515 double precision p1,esoll
0516 dimension p1(5),p0(5,mamx+mamx)
0517 call utpri('utrsph',ish,ishini,4)
0518
0519 errlim=0.0001
0520
0521 iret=0
0522 nptlpt=iabs(maproj)+iabs(matarg)
0523 if(nptl.le.nptlpt) goto 9999
0524
0525 esoll=0.d0
0526 p1(1)=0.d0
0527 p1(2)=0.d0
0528 p1(3)=0.d0
0529 p1(4)=0.d0
0530 do i=nptlpt+1,nptl
0531 if((istptl(i).le.11
0532 $ .and.(iorptl(i).ge.1.and.istptl(iorptl(i)).eq.41))
0533 $ .or.istptl(i).eq.20.or.istptl(i).eq.21)then
0534 do j=1,2
0535 p1(j)=p1(j)+dble(pptl(j,i))
0536 enddo
0537 endif
0538 enddo
0539 do i=1,nptlpt
0540 do j=1,5
0541 p0(j,i)=pptl(j,i)
0542 enddo
0543 do j=3,4
0544 p1(j)=p1(j)+dble(pptl(j,i))
0545 enddo
0546 enddo
0547 p1(5)=dsqrt((p1(4)+p1(3))*(p1(4)-p1(3))-p1(2)**2.d0-p1(1)**2.d0)
0548 esoll=p1(5)
0549 if(ish.ge.4) write (ifch,'(a,5g13.6)') 'boost-vector',p1
0550
0551 c trafo
0552 c -----
0553 do i=1,nptl
0554 if((istptl(i).le.11
0555 $ .and.(iorptl(i).ge.1.and.istptl(iorptl(i)).eq.41))
0556 $ .or.istptl(i).eq.20.or.istptl(i).eq.21
0557 $ .or.(istptl(i).eq.0.and.i.le.nptlpt))then
0558 call utlob4(1,p1(1),p1(2),p1(3),p1(4),p1(5)
0559 $ ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i))
0560 endif
0561 enddo
0562
0563
0564 if(ish.ge.5)write(ifch,'(a)')'------------------'
0565
0566 c rescale momenta in rest frame
0567 c -----------------------------
0568
0569 scal=1.
0570 diff=0.
0571 do ipass=1,1000
0572 sum=0.
0573 sum3=0.
0574 ndif=0
0575 do j=1,nptl
0576 if((istptl(j).le.11
0577 $ .and.(iorptl(j).ge.1.and.istptl(iorptl(j)).eq.41))
0578 $ .or.istptl(j).eq.20.or.istptl(j).eq.21
0579 $ .or.(istptl(j).eq.0.and.j.le.nptlpt))then
0580 if(j.gt.nptlpt)then
0581 ndif=ndif+1
0582 pptl(3,j)=scal*(pptl(3,j)-diff)
0583 pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2
0584 * +pptl(5,j)**2)
0585 endif
0586 sum=sum+pptl(4,j)
0587 sum3=sum3+pptl(3,j)
0588 endif
0589 enddo
0590
0591 diff=sum3/real(ndif)
0592 scal=real(esoll)/sum
0593 if(ish.ge.6)write(ifch,*)'ipass,scal,diff,e,esoll,pz,ndif:'
0594 $ ,ipass,scal,diff,sum,esoll,sum3,ndif
0595 if(abs(scal-1.).le.errlim.and.abs(diff).lt.10.*errlim) goto300
0596 enddo
0597 if(ish.ge.1)then
0598 call utmsg('hresph')
0599 write(ifch,*)'***** scal=',scal,diff
0600 call utmsgf
0601 endif
0602
0603
0604 c trafo
0605 c -----
0606 300 continue
0607 c do i=nptlpt+1,nptl
0608 do i=1,nptl
0609 if((istptl(i).le.11
0610 $ .and.(iorptl(i).ge.1.and.istptl(iorptl(i)).eq.41))
0611 $ .or.istptl(i).eq.20.or.istptl(i).eq.21
0612 $ .or.(istptl(i).eq.0.and.i.le.nptlpt))then
0613 call utlob4(-1,p1(1),p1(2),p1(3),p1(4),p1(5)
0614 $ ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i))
0615 endif
0616 if(i.le.nptlpt)then
0617 do j=1,5
0618 pptl(j,i)=p0(j,i)
0619 enddo
0620 endif
0621 enddo
0622
0623 9999 call utprix('utrsph',ish,ishini,4)
0624
0625 end
0626
0627 cc-----------------------------------------------------------------------
0628 c double precision function dddlog(xxx)
0629 cc-----------------------------------------------------------------------
0630 c double precision xxx
0631 c dddlog=-1d50
0632 c if(xxx.gt.0d0)dddlog=dlog(xxx)
0633 c end
0634 c
0635 ccc-----------------------------------------------------------------------
0636 c subroutine randfl(jc,iqa0,iflav,ic,isame)
0637 cc-----------------------------------------------------------------------
0638 cc returns random flavour ic(2) (iqa0=1:quark,2:antiquark,11:diquark)
0639 cc-----------------------------------------------------------------------
0640 c include 'epos.inc'
0641 c real probab(nflav),probsu(nflav+1)
0642 c integer jc(nflav,2),jc0(nflav,2),ic(2)
0643 c if(ish.ge.6)then
0644 c write(ifch,*)('-',i=1,10)
0645 c *,' entry sr randfl ',('-',i=1,30)
0646 c write(ifch,*)'iqa0:',iqa0
0647 c write(ifch,*)'jc:'
0648 c write(ifch,*)jc
0649 c endif
0650 c iflav=0
0651 c ic(1)=0
0652 c ic(2)=0
0653 c do 10 n=1,nflav
0654 c do 10 i=1,2
0655 c10 jc0(n,i)=0
0656 c iqa1=iqa0*10
0657 c9999 iqa1=iqa1/10
0658 c if(iqa1.eq.0)goto9998
0659 c iqa=mod(iqa1,10)
0660 c su=0
0661 c do 20 i=1,nflav
0662 c probab(i)=jc(i,iqa)-jc0(i,iqa)
0663 c if(isame.eq.1)probab(i)=probab(i)*(jc(i,3-iqa)-jc0(i,3-iqa))
0664 c20 su=su+probab(i)
0665 c if(su.lt..5)then
0666 c iflav=0
0667 c ic(1)=0
0668 c ic(2)=0
0669 c goto9998
0670 c endif
0671 c probsu(1)=0.
0672 c do 30 i=1,nflav
0673 c probsu(i+1)=probsu(i)+probab(i)/su
0674 c if(probsu(i+1)-probsu(i).lt.1e-5)probsu(i+1)=probsu(i)
0675 c30 continue
0676 c r=rangen()*probsu(nflav+1)
0677 c do 50 i=1,nflav
0678 c if(probsu(i).le.r.and.r.lt.probsu(i+1))iflav=i
0679 c50 continue
0680 c jc0(iflav,iqa)=jc0(iflav,iqa)+1
0681 c if(isame.eq.1)jc0(iflav,3-iqa)=jc0(iflav,3-iqa)+1
0682 c call idenco(jc0,ic,ireten)
0683 c if(ireten.eq.1)call utstop('randfl: idenco ret code = 1&')
0684 c if(ish.ge.6)then
0685 c write(ifch,*)'probab:'
0686 c write(ifch,*)probab
0687 c write(ifch,*)'probsu:'
0688 c write(ifch,*)probsu
0689 c write(ifch,*)'ran#:',r,' flav:',iflav
0690 c endif
0691 c goto9999
0692 c9998 continue
0693 c if(ish.ge.6)write(ifch,*)('-',i=1,30)
0694 c *,' exit sr randfl ',('-',i=1,10)
0695 c return
0696 c end
0697 c
0698 c
0699 cc-----------------------------------------------------------------------
0700 c subroutine ranhvy(x,eps)
0701 cc-----------------------------------------------------------------------
0702 cc generates x for heavy particle fragmentation according to
0703 cc the peterson form
0704 cc d(x)=1/(x*(1-1/x-eps/(1-x))**2)
0705 cc =d0(x)*d1(x)*d2(x)
0706 cc d0(x)=(1-x)**2/((1-x)**2+eps)**2
0707 cc d1(x)=x
0708 cc d2(x)=(((1-x)**2+eps)/((1-x)**2+eps*x))**2
0709 cc using x=1-y**pow
0710 cc generates flat in x if eps>1.
0711 cc-----------------------------------------------------------------------
0712 c data aln4/1.3863/
0713 c if(eps.lt.1.) then
0714 c pow=alog((3.+eps)/eps)/aln4
0715 c ymx=(eps*(3.*pow-1.)/(pow+1.))**(.5/pow)
0716 c zmx=1-ymx**pow
0717 c d0mx=(1-zmx)**2/((1.-zmx)**2+eps)**2*pow*ymx**(pow-1.)
0718 c d2mx=2./(2.-sqrt(eps))
0719 c else
0720 c pow=1.
0721 c zmx=0.
0722 c d0mx=(1.-zmx)**2/((1.-zmx)**2+eps)**2
0723 c d2mx=1.+eps
0724 c endif
0725 cc
0726 cc generate z according to (1-z)**2/((1-z)**2+eps*z)**2
0727 c1 continue
0728 c y=rangen()
0729 c z=1.-y**pow
0730 cc
0731 c d0z=(1.-z)**2/((1.-z)**2+eps)**2*pow*y**(pow-1.)
0732 c if(d0z.lt.rangen()*d0mx) goto1
0733 cc
0734 cc check remaining factors
0735 c d1=z
0736 c d2=(((1.-z)**2+eps)/((1.-z)**2+eps*z))**2
0737 c if(d1*d2.lt.rangen()*d2mx) goto1
0738 cc
0739 cc good x
0740 c x=z
0741 c return
0742 c end
0743 c
0744 *-- Author : D. HECK IK FZK KARLSRUHE 27/04/1994
0745 C=======================================================================
0746
0747 SUBROUTINE EPOVAPOR( MAPRO,INEW,JFIN,ITYP,PFRX,PFRY,PFRZ )
0748
0749 C-----------------------------------------------------------------------
0750 C (E)VAPOR(ATION OF NUCLEONS AND ALPHA PARTICLES FROM FRAGMENT)
0751 C
0752 C TREATES THE REMAINING UNFRAGMENTED NUCLEUS
0753 C EVAPORATION FOLLOWING CAMPI APPROXIMATION.
0754 C SEE: X. CAMPI AND J. HUEFNER, PHYS.REV. C24 (1981) 2199
0755 C AND J.J. GAIMARD, THESE UNIVERSITE PARIS 7, (1990)
0756 C THIS SUBROUTINE IS CALLED FROM SDPM, DPMJST, NSTORE, AND VSTORE.
0757 C ARGUMENTS INPUT:
0758 C MAPRO = NUMBER OF NUCLEONS OF PROJECTILE
0759 C INEW = PARTICLE TYPE OF SPECTATOR FRAGMENT
0760 C ARGUMENTS OUTPUT:
0761 C JFIN = NUMBER OF FRAGMENTS
0762 C ITYP(1:JFIN) = NATURE (PARTICLE CODE) OF FRAGMENTS
0763 C PFRX(1:JFIN) = TRANSVERSE MOMENTUM OF FRAGMENTS IN X-DIRECTION
0764 C PFRY(1:JFIN) = TRANSVERSE MOMENTUM OF FRAGMENTS IN Y-DIRECTION
0765 C PFRZ(1:JFIN) = LONGITUDINAL MOMENTUM OF FRAGMENTS IN Z-DIRECTION
0766 C
0767 C FROM CORSIKA AND ADAPTED BY T. PIEROG TO INCLUDE LONG MOMENTUM AND
0768 C MORE REALISTIC FRAGMENTS
0769 C-----------------------------------------------------------------------
0770
0771 IMPLICIT NONE
0772 include 'epos.inc'
0773 common/eporansto2/irndmseq
0774 integer irndmseq
0775
0776 DOUBLE PRECISION PFR(mamxx),PFRX(mamxx),PFRY(mamxx),PFRZ(mamxx)
0777 * ,RD(2*mamxx),SPFRY,SPFRZ,drangen
0778 DOUBLE PRECISION AFIN,AGLH,APRF,BGLH,EEX,PHIFR,RANNORM,SPFRX
0779 INTEGER ITYP(mamxx),IARM,INEW,ITYPRM,INRM,IS,IZRM,JC,lseq
0780 * ,JFIN,K,L,LS,MAPRO,MF,NFIN,NINTA,NNUC,NPRF,NNSTEP
0781 SAVE
0782 EXTERNAL RANNORM
0783 C-----------------------------------------------------------------------
0784
0785 IF(ish.ge.7)WRITE(ifch,*)'EPOVAPOR : MAPRO,INEW=',MAPRO,INEW
0786
0787 lseq = irndmseq
0788 ITYPRM = INEW
0789 NPRF = INEW/100
0790 NINTA = MAPRO - NPRF
0791 IF ( NINTA .EQ. 0 ) THEN
0792 C NO NUCLEON HAS INTERACTED
0793 JFIN = 1
0794 PFRX(1) = 0.D0
0795 PFRY(1) = 0.D0
0796 PFRZ(1) = 0.D0
0797 ITYP(1) = -ITYPRM
0798 IF(ish.ge.7)WRITE(ifch,*) 'EPOVAPOR : JFIN,NINTA=',JFIN,NINTA
0799 RETURN
0800 ENDIF
0801
0802 C EXCITATION ENERGY EEX OF PREFRAGMENT
0803 C SEE: J.J. GAIMARD, THESE UNIVERSITE PARIS 7, (1990), CHPT. 4.2
0804 EEX = 0.D0
0805 CALL RMMARD( RD,2*NINTA,lseq )
0806 DO L = 1, NINTA
0807 IF ( RD(NINTA+L) .LT. RD(L) ) RD(L) = 1.D0 - RD(L)
0808 EEX = EEX + RD(L)
0809 ENDDO
0810 C DEPTH OF WOODS-SAXON POTENTIAL TO FERMI SURFACE IS 0.040 GEV
0811 IF(ish.ge.7)WRITE(ifch,*)'EPOVAPOR : EEX=',SNGL(EEX*0.04D0),
0812 & ' GEV'
0813 C EVAPORATION: EACH EVAPORATION STEP NEEDS ABOUT 0.020 GEV, THEREFORE
0814 C NNSTEP IS EEX * 0.04/0.02 = EEX * 2.
0815 NNSTEP = INT( EEX*2.D0 )
0816
0817 IF ( NNSTEP .LE. 0 ) THEN
0818 C EXCITATION ENERGY TOO SMALL, NO EVAPORATION
0819 JFIN = 1
0820 PFRX(1) = 0.D0
0821 PFRY(1) = 0.D0
0822 PFRZ(1) = 0.D0
0823 ITYP(1) = -ITYPRM
0824 IF(ish.ge.7)WRITE(ifch,*) 'EPOVAPOR : JFIN,EEX=',JFIN,SNGL(EEX)
0825 RETURN
0826 ENDIF
0827
0828 C AFIN IS ATOMIC NUMBER OF FINAL NUCLEUS
0829 APRF = DBLE(NPRF)
0830 AFIN = APRF - 1.6D0 * DBLE(NNSTEP)
0831 NFIN = MAX( 0, INT( AFIN+0.5D0 ) )
0832 C CORRESPONDS TO DEFINITION; FRAGMENTATION-EVAPORATION
0833 C CONVOLUTION EMU07 /MODEL ABRASION EVAPORATION (JNC FZK APRIL 94)
0834 C NNUC IS NUMBER OF EVAPORATING NUCLEONS
0835 NNUC = NPRF - NFIN
0836 IF(ish.ge.7)WRITE(ifch,*) 'EPOVAPOR : NFIN,NNUC=',NFIN,NNUC
0837 JC = 0
0838
0839 IF ( NNUC .LE. 0 ) THEN
0840 C NO EVAPORATION
0841 JFIN = 1
0842 PFRX(1) = 0.D0
0843 PFRY(1) = 0.D0
0844 PFRZ(1) = 0.D0
0845 ITYP(1) = -ITYPRM
0846 RETURN
0847
0848 ELSEIF ( NNUC .GE. 4 ) THEN
0849 C EVAPORATION WITH FORMATION OF ALPHA PARTICLES POSSIBLE
0850 C IARM, IZRM, INRM ARE NUMBER OF NUCLEONS, PROTONS, NEUTRONS OF
0851 C REMAINDER
0852 DO LS = 1, NNSTEP
0853 IARM = ITYPRM/100
0854 IF ( IARM .LE. 0 ) GOTO 100
0855 IZRM = MOD(ITYPRM,100)
0856 INRM = IARM - IZRM
0857 JC = JC + 1
0858 CALL RMMARD( RD,2,lseq )
0859 IF ( RD(1) .LT. 0.2D0 .AND. IZRM .GE. 2
0860 * .AND. INRM .GE. 2 ) THEN
0861 ITYP(JC) = -402 !alpha
0862 NNUC = NNUC - 4
0863 ITYPRM = ITYPRM - 402
0864 ELSE
0865 IF ( IZRM .EQ. 1 .AND. INRM .GT. IZRM ) THEN
0866 ITYP(JC) = 1220
0867 ITYPRM = ITYPRM - 100
0868 ELSEIF ( INRM .EQ. 1 .AND. IZRM .GT. INRM ) THEN
0869 ITYP(JC) = 1120
0870 ITYPRM = ITYPRM - 101
0871 ELSEIF(RD(2)*(IZRM+INRM).LT.IZRM.AND.IZRM.GE.INRM)THEN
0872 ITYP(JC) = 1120
0873 ITYPRM = ITYPRM - 101
0874 ELSE
0875 ITYP(JC) = 1220
0876 ITYPRM = ITYPRM - 100
0877 ENDIF
0878 NNUC = NNUC - 1
0879 ENDIF
0880 IF ( NNUC .LE. 0 ) GOTO 50
0881 ENDDO
0882 ENDIF
0883
0884 IF ( NNUC .LT. 4 ) THEN
0885 C EVAPORATION WITHOUT FORMATION OF ALPHA PARTICLES
0886 CALL RMMARD( RD,NNUC,lseq )
0887 DO IS = 1, NNUC
0888 IARM = ITYPRM/100
0889 IF ( IARM .LE. 0 ) GOTO 100
0890 IZRM = MOD(ITYPRM,100)
0891 INRM = IARM - IZRM
0892 JC = JC + 1
0893 IF ( IZRM .EQ. 1 .AND. INRM .GT. IZRM ) THEN
0894 ITYP(JC) = 1220
0895 ITYPRM = ITYPRM - 100
0896 ELSEIF ( INRM .EQ. 1 .AND. IZRM .GT. IZRM ) THEN
0897 ITYP(JC) = 1120
0898 ITYPRM = ITYPRM - 101
0899 ELSEIF ( RD(IS)*IARM .LT. IZRM .AND. IZRM .GE. INRM ) THEN
0900 ITYP(JC) = 1120
0901 ITYPRM = ITYPRM - 101
0902 ELSE
0903 ITYP(JC) = 1220
0904 ITYPRM = ITYPRM - 100
0905 ENDIF
0906 ENDDO
0907 ENDIF
0908
0909 50 CONTINUE
0910 IARM = ITYPRM/100
0911 IF ( IARM .LE. 0 ) GOTO 100
0912 IZRM = MOD(ITYPRM,100)
0913 INRM = IARM - IZRM
0914 JC = JC + 1
0915 IF ( IARM .EQ. 5 ) THEN !EXCLUDED
0916 IF ( IZRM .GE. INRM ) THEN
0917 ITYP(JC) = 1120
0918 ITYPRM = ITYPRM - 101
0919 ELSE
0920 ITYP(JC) = 1220
0921 ITYPRM = ITYPRM - 100
0922 ENDIF
0923 JC = JC + 1
0924 ITYP(JC) = -ITYPRM
0925 ELSEIF ( ITYPRM .GT. 200 ) THEN
0926 ITYP(JC) = -ITYPRM
0927 ELSEIF ( ITYPRM .EQ. 200 ) THEN
0928 ITYP(JC) = 1220
0929 JC = JC + 1
0930 ITYP(JC) = 1220
0931 ELSEIF ( ITYPRM .EQ. 101 ) THEN
0932 ITYP(JC) = 1120
0933 ELSEIF ( ITYPRM .EQ. 100 ) THEN
0934 ITYP(JC) = 1220
0935 ELSE
0936 JC = JC - 1
0937 IF ( ITYPRM .NE. 0 ) WRITE(*,*)
0938 * 'EPOVAPOR : ILLEGAL PARTICLE ITYPRM =',ITYPRM
0939 ENDIF
0940
0941 100 CONTINUE
0942 JFIN = JC
0943 IF(ish.ge.7)WRITE(ifch,*)
0944 * 'EPOVAPOR : NO ITYP PFR PFL'
0945 IF ( infragm .EQ. 2 ) THEN
0946 C EVAPORATION WITH PT AFTER PARAMETRIZED JACEE DATA
0947 DO MF = 1, JFIN
0948 IF(ITYP(MF).LT.0)THEN
0949 IARM=-ITYP(MF)/100
0950 ELSE
0951 IARM=1
0952 ENDIF
0953 PFR(MF) = RANNORM(0.088D0,0.044D0)
0954 PFRZ(MF)= (2*int(drangen(PFR(MF))+0.5d0)-1)
0955 & *RANNORM(0.300D0/DBLE(IARM),0.100D0/SQRT(DBLE(IARM))) !Fermi motion about 300 MeV
0956 IF(ish.ge.7)WRITE(ifch,*) MF,ITYP(MF),SNGL(PFR(MF))
0957 & ,SNGL(PFRZ(MF))
0958 ENDDO
0959 ELSEIF ( infragm .EQ. 3 ) THEN
0960 C EVAPORATION WITH PT AFTER GOLDHABER''S MODEL (PHYS.LETT.53B(1974)306)
0961 DO MF = 1, JFIN
0962 K = MAX( 1, -ITYP(MF)/100 )
0963 BGLH = K * (MAPRO - K) / DBLE(MAPRO-1)
0964 C THE VALUE 0.103 [GEV] IS SIGMA(0)=P(FERMI)/SQRT(5.)
0965 * AGLH = 0.103D0 * SQRT( BGLH )
0966 C THE VALUE 0.090 [GEV] IS EXPERIMENTALLY DETERMINED SIGMA(0)
0967 AGLH = 0.090D0 * SQRT( BGLH )
0968 PFR(MF) = RANNORM(0.D0,AGLH)
0969 PFRZ(MF)= RANNORM(0.000D0,0.500D0) !from pAg at 100 GeV
0970 IF(ish.ge.7)WRITE(ifch,*) MF,ITYP(MF),SNGL(PFR(MF))
0971 & ,SNGL(PFRZ(MF))
0972 ENDDO
0973 ELSE
0974 C EVAPORATION WITHOUT TRANSVERSE MOMENTUM
0975 DO MF = 1, JFIN
0976 PFR(MF) = 0.D0
0977 PFRZ(MF)= 0.D0
0978 IF(ish.ge.7)WRITE(ifch,*) MF,ITYP(MF),SNGL(PFR(MF))
0979 & ,SNGL(PFRZ(MF))
0980 ENDDO
0981 ENDIF
0982 C CALCULATE RESIDUAL TRANSVERSE MOMENTUM
0983 SPFRX = 0.D0
0984 SPFRY = 0.D0
0985 SPFRZ = 0.D0
0986 CALL RMMARD( RD,JFIN,lseq )
0987 DO MF = 1, JFIN
0988 PHIFR = PI * RD(MF)
0989 PFRX(MF) = PFR(MF) * COS( PHIFR )
0990 PFRY(MF) = PFR(MF) * SIN( PHIFR )
0991 SPFRY = SPFRY + PFRY(MF)
0992 SPFRX = SPFRX + PFRX(MF)
0993 SPFRZ = SPFRZ + PFRZ(MF)
0994 ENDDO
0995 C CORRECT ALL TRANSVERSE MOMENTA FOR MOMENTUM CONSERVATION
0996 SPFRX = SPFRX / JFIN
0997 SPFRY = SPFRY / JFIN
0998 SPFRZ = SPFRZ / JFIN
0999 DO MF = 1, JFIN
1000 PFRX(MF) = PFRX(MF) - SPFRX
1001 PFRY(MF) = PFRY(MF) - SPFRY
1002 PFRZ(MF) = PFRZ(MF) - SPFRZ
1003 ENDDO
1004
1005 IF(ish.ge.7)WRITE(ifch,*) 'EPOVAPOR : NINTA,JFIN=',NINTA,JFIN
1006
1007 RETURN
1008 END
1009
1010 *-- Author : The CORSIKA development group 21/04/1994
1011 C=======================================================================
1012
1013 DOUBLE PRECISION FUNCTION RANNORM( A,B )
1014
1015 C-----------------------------------------------------------------------
1016 C RAN(DOM NUMBER) NOR(MALLY DISTRIBUTED)
1017 C
1018 C GENERATES NORMAL DISTRIBUTED RANDOM NUMBER
1019 C DELIVERS 2 UNCORRELATED RANDOM NUMBERS,
1020 C THEREFORE RANDOM CALLS ARE ONLY NECESSARY EVERY SECOND TIME.
1021 C REFERENCE : NUMERICAL RECIPES, W.H. PRESS ET AL.,
1022 C CAMBRIDGE UNIVERSITY PRESS, 1992 ISBN 0 521 43064 X
1023 C ARGUMENTS:
1024 C A = MEAN VALUE
1025 C B = STANDARD DEVIATION
1026 C
1027 C FROM CORSIKA
1028 C-----------------------------------------------------------------------
1029
1030 IMPLICIT NONE
1031 double precision facrdm,u1rdm,u2rdm,drangen
1032 logical knordm
1033 data knordm/.true./
1034
1035 DOUBLE PRECISION A,B,RR
1036 SAVE facrdm,u1rdm,u2rdm,knordm
1037 C-----------------------------------------------------------------------
1038
1039 IF ( KNORdm ) THEN
1040 1 CONTINUE
1041 U1rdm = 2.D0*drangen(a) - 1.D0
1042 U2rdm = 2.D0*drangen(b) - 1.D0
1043 RR = U1rdm**2 + U2rdm**2
1044 IF ( RR .GE. 1.D0 .OR. RR .EQ. 0.D0 ) GOTO 1
1045 FACrdm = SQRT( (-2.D0) * LOG(RR) / RR )
1046
1047 RANNORM = FACrdm * U1rdm * B + A
1048 KNORdm = .FALSE.
1049 ELSE
1050 RANNORM = FACrdm * U2rdm * B + A
1051 KNORdm = .TRUE.
1052 ENDIF
1053
1054 RETURN
1055 END
1056
1057 c-----------------------------------------------------------------------
1058 function ransig()
1059 c-----------------------------------------------------------------------
1060 c returns randomly +1 or -1
1061 c-----------------------------------------------------------------------
1062 ransig=1
1063 if(rangen().gt.0.5)ransig=-1
1064 return
1065 end
1066
1067 cc-----------------------------------------------------------------------
1068 c function ranxq(n,x,q,xmin)
1069 cc-----------------------------------------------------------------------
1070 cc returns random number according to x(i) q(i) with x>=xmin
1071 cc-----------------------------------------------------------------------
1072 c include 'epos.inc'
1073 c real x(n),q(n)
1074 c imin=1
1075 c if(xmin.eq.0.)goto3
1076 c i1=1
1077 c i2=n
1078 c1 i=i1+(i2-i1)/2
1079 c if(x(i).lt.xmin)then
1080 c i1=i
1081 c elseif(x(i).gt.xmin)then
1082 c i2=i
1083 c else
1084 c imin=i
1085 c goto3
1086 c endif
1087 c if(i2-i1.gt.1)goto1
1088 c imin=i2
1089 c3 continue
1090 c if(q(imin).gt.q(n)*.9999)then
1091 c ranxq=xmin
1092 c goto4
1093 c endif
1094 c qran=q(imin)+rangen()*(q(n)-q(imin))
1095 c ranxq=utinvt(n,x,q,qran)
1096 c4 continue
1097 c
1098 c if(ranxq.lt.xmin)then
1099 c call utmsg('ranxq ')
1100 c write(ifch,*)'***** ranxq=',ranxq,' < xmin=',xmin
1101 c write(ifch,*)'q(imin) q q(n):',q(imin),qran,q(n)
1102 c write(ifch,*)'x(imin) x x(n):',x(imin),ranxq,x(n)
1103 c call utmsgf
1104 c ranxq=xmin
1105 c endif
1106 c
1107 c return
1108 c end
1109 c
1110 cc ***** end r-routines
1111 cc ***** beg s-routines
1112 c
1113 cc-----------------------------------------------------------------------
1114 c function sbet(z,w)
1115 cc-----------------------------------------------------------------------
1116 c sbet=utgam1(z)*utgam1(w)/utgam1(z+w)
1117 c return
1118 c end
1119 c
1120 cc-----------------------------------------------------------------------
1121 c function smass(a,y,z)
1122 cc-----------------------------------------------------------------------
1123 cc returns droplet mass (in gev) (per droplet, not (!) per nucleon)
1124 cc according to berger/jaffe mass formula, prc35(1987)213 eq.2.31,
1125 cc see also c. dover, BNL-46322, intersections-meeting, tucson, 91.
1126 cc a: massnr, y: hypercharge, z: charge,
1127 cc-----------------------------------------------------------------------
1128 c common/cmass/thet,epsi,as,ac,dy,dz,ym,cz,zm,sigma,rzero
1129 c ymin=ym*a
1130 c zmin=cz/(dz/a+zm/a**(1./3.))
1131 c smass=epsi*a+as*a**(2./3.)+(ac/a**(1./3.)+dz/a/2.)*(z-zmin)**2
1132 c *+dy/a/2.*(y-ymin)**2
1133 c return
1134 c end
1135 c
1136 cc-----------------------------------------------------------------------
1137 c subroutine smassi(theta)
1138 cc-----------------------------------------------------------------------
1139 cc initialization for smass.
1140 cc calculates parameters for berger/jaffe mass formula
1141 cc (prc35(1987)213 eq.2.31, see also c. dover, BNL-46322).
1142 cc theta: parameter that determines all parameters in mass formula.
1143 cc-----------------------------------------------------------------------
1144 c common/cmass/thet,epsi,as,ac,dy,dz,ym,cz,zm,sigma,rzero
1145 c thet=theta
1146 c
1147 c astr=.150
1148 c pi=3.14159
1149 c alp=1./137.
1150 c
1151 c co=cos(theta)
1152 c si=sin(theta)
1153 c bet=(1+co**3)/2.
1154 c rzero=si/astr/( 2./3./pi*(1+co**3) )**(1./3.)
1155 cctp060829 cs=astr/si
1156 c cz=-astr/si*( ( .5*(1+co**3) )**(1./3.)-1 )
1157 c sigma=6./8./pi*(astr/si)**3*(co**2/6.-si**2*(1-si)/3.-
1158 c *1./3./pi*(pi/2.-theta-sin(2*theta)+si**3*alog((1+co)/si)))
1159 c
1160 c epsi=astr*((.5*(1+co**3))**(1./3.)+2)/si
1161 c as=4*pi*sigma*rzero**2
1162 c ac=3./5.*alp/rzero
1163 c dz=astr/si*bet**(1./3.)*co**2*
1164 c *(co**4*(1+bet**(2./3.))+(1+bet)**2)/
1165 c *( (2*co**2+bet**(1./3.))*(co**4*(1+bet**(2./3.))+(1+bet)**2)-
1166 c *(co**4+bet**(1./3.)*(1+bet))*((2*bet**(2./3.)-1)*co**2+1+bet) )
1167 c dy=astr/6.*(1+co**3)**3/si*
1168 c *( 1+(1+co)/(4*(1+co**3))**(2./3.) )/
1169 c *(co**6+co+co*(.5*(1+co**3))**(4./3.))
1170 c zm=6*alp/(5*rzero)
1171 c ym=(1-co**3)/(1+co**3)
1172 c
1173 c return
1174 c end
1175 c
1176 cc-----------------------------------------------------------------------
1177 c subroutine smassp
1178 cc-----------------------------------------------------------------------
1179 cc prints smass.
1180 cc-----------------------------------------------------------------------
1181 c include 'epos.inc'
1182 c common/cmass/thet,epsi,as,ac,dy,dz,ym,cz,zm,sigma,rzero
1183 c real eng(14),ymi(14),zmi(14)
1184 c pi=3.14159
1185 c write(ifch,*)'parameters of mass formula:'
1186 c write(ifch,*)'theta=',thet,' epsi=',epsi
1187 c write(ifch,*)'as=',as,' ac=',ac
1188 c write(ifch,*)'dy=',dy,' dz=',dz
1189 c write(ifch,*)'ym=',ym
1190 c write(ifch,*)'cz dz zm=',cz,dz,zm
1191 c write(ifch,*)'sigma**1/3=',sigma**(1./3.),' rzero=',rzero
1192 c write(ifch,*)'mass:'
1193 c write(ifch,5000)(j,j=1,14)
1194 c5000 format(5x,'a:',14i5)
1195 c do 4 j=1,14
1196 c a=j
1197 c ymi(j)=ym*a
1198 c4 zmi(j)=cz/(dz/a+zm/a**(1./3.))
1199 c write(ifch,5002)(ymi(j),j=1,14)
1200 c5002 format(1x,'ymin: ',14f5.2)
1201 c write(ifch,5003)(zmi(j),j=1,14)
1202 c5003 format(1x,'zmin: ',14f5.2)
1203 c do 2 i=1,15
1204 c ns=11-i
1205 c do 3 j=1,14
1206 c a=j
1207 c y=a-ns
1208 c z=0.
1209 c3 eng(j)=smass(a,y,z)/a
1210 c write(ifch,5001)ns,(eng(j),j=1,14)
1211 c5001 format(1x,'s=',i2,2x,14f5.2)
1212 c2 continue
1213 c write(ifch,*)'mass-mass(free):'
1214 c write(ifch,5000)(j,j=1,14)
1215 c do 5 i=1,15
1216 c ns=11-i
1217 c do 6 j=1,14
1218 c a=j
1219 c y=a-ns
1220 c z=0.
1221 c call smassu(a,y,z,ku,kd,ks,kc)
1222 c6 eng(j)=(smass(a,y,z)-utamnu(ku,kd,ks,kc,0,0,3))/a
1223 c write(ifch,5001)ns,(eng(j),j=1,14)
1224 c5 continue
1225 c
1226 c stop
1227 c end
1228 c
1229 cc-----------------------------------------------------------------------
1230 c subroutine smasst(kux,kdx,ksx,kcx,a,y,z)
1231 cc-----------------------------------------------------------------------
1232 cc input: kux,kdx,ksx,kcx = net quark numbers (for u,d,s,c quarks).
1233 cc output: massnr a, hypercharge y and charge z.
1234 cc-----------------------------------------------------------------------
1235 c sg=1
1236 c if(kux+kdx+ksx+kcx.lt.0.)sg=-1
1237 c ku=sg*kux
1238 c kd=sg*kdx
1239 c ks=sg*ksx
1240 c kc=sg*kcx
1241 c k=ku+kd+ks+kc
1242 c if(mod(k,3).ne.0)stop'noninteger baryon number'
1243 c a=k/3
1244 c y=a-ks
1245 c nz=2*ku-kd-ks+2*kc
1246 c if(mod(nz,3).ne.0)stop'noninteger charge'
1247 c z=nz/3
1248 c return
1249 c end
1250 c
1251 cc-----------------------------------------------------------------------
1252 c subroutine smassu(ax,yx,zx,ku,kd,ks,kc)
1253 cc-----------------------------------------------------------------------
1254 cc input: massnr ax, hypercharge yx and charge zx.
1255 cc output: ku,kd,ks,kc = net quark numbers (for u,d,s,c quarks).
1256 cc-----------------------------------------------------------------------
1257 c sg=1
1258 c if(ax.lt.0.)sg=-1
1259 c a=sg*ax
1260 c y=sg*yx
1261 c z=sg*zx
1262 c ku=nint(a+z)
1263 c kd=nint(a-z+y)
1264 c ks=nint(a-y)
1265 c kc=0
1266 c return
1267 c end
1268 c
1269 cc-----------------------------------------------------------------------
1270 c function spoc(a,b,c,d,x)
1271 cc-----------------------------------------------------------------------
1272 cc power fctn with cutoff
1273 cc-----------------------------------------------------------------------
1274 c spoc=0
1275 c if(a.eq.0..and.b.eq.0.)return
1276 c spoc =a+b*x**c
1277 c spoc0=a+b*d**c
1278 c spoc=amin1(spoc,spoc0)
1279 c spoc=amax1(0.,spoc)
1280 c return
1281 c end
1282 c
1283 c-----------------------------------------------------------------------
1284 function utacos(x)
1285 c-----------------------------------------------------------------------
1286 c returns acos(x) for -1 <= x <= 1 , acos(+-1) else
1287 c-----------------------------------------------------------------------
1288 include 'epos.inc'
1289 argum=x
1290 if(x.lt.-1.)then
1291 if(ish.ge.1)then
1292 call utmsg('utacos')
1293 write(ifch,*)'***** argum = ',argum,' set -1'
1294 call utmsgf
1295 endif
1296 argum=-1.
1297 elseif(x.gt.1.)then
1298 if(ish.ge.1)then
1299 call utmsg('utacos')
1300 write(ifch,*)'***** argum = ',argum,' set 1'
1301 call utmsgf
1302 endif
1303 argum=1.
1304 endif
1305 utacos=acos(argum)
1306 return
1307 end
1308
1309 c----------------------------------------------------------------------
1310 function utamnu(keux,kedx,kesx,kecx,kebx,ketx,modus)
1311 c----------------------------------------------------------------------
1312 c returns min mass of droplet with given u,d,s,c content
1313 c keux: net u quark number
1314 c kedx: net d quark number
1315 c kesx: net s quark number
1316 c kecx: net c quark number
1317 c kebx: net b quark number
1318 c ketx: net t quark number
1319 c modus: 4=two lowest multiplets; 5=lowest multiplet
1320 c----------------------------------------------------------------------
1321 common/files/ifop,ifmt,ifch,ifcx,ifhi,ifdt,ifcp,ifdr
1322 common/csjcga/amnull,asuha(7)
1323 common/drop4/asuhax(7),asuhay(7)
1324
1325 if(modus.lt.4.or.modus.gt.5)stop'UTAMNU: not supported'
1326 c 1 format(' flavours:',6i5 )
1327 c 100 format(' flavours+mass:',6i5,f8.2 )
1328 c write(ifch,1)keux,kedx,kesx,kecx,kebx,ketx
1329
1330 amnull=0.
1331
1332 do i=1,7
1333 if(modus.eq.4)asuha(i)=asuhax(i) !two lowest multiplets
1334 if(modus.eq.5)asuha(i)=asuhay(i) !lowest multiplet
1335 enddo
1336
1337 ke=iabs(keux+kedx+kesx+kecx+kebx+ketx)
1338
1339 if(keux+kedx+kesx+kecx+kebx+ketx.ge.0)then
1340 keu=keux
1341 ked=kedx
1342 kes=kesx
1343 kec=kecx
1344 keb=kebx
1345 ket=ketx
1346 else
1347 keu=-keux
1348 ked=-kedx
1349 kes=-kesx
1350 kec=-kecx
1351 keb=-kebx
1352 ket=-ketx
1353 endif
1354
1355 c write(ifch,*)keu,ked,kes,kec,keb,ket
1356
1357 c removing top mesons to remove t quarks or antiquarks
1358 if(ket.ne.0)then
1359 12 continue
1360 ii=sign(1,ket)
1361 ket=ket-ii
1362 if(ii*keu.le.ii*ked)then
1363 keu=keu+ii
1364 else
1365 ked=ked+ii
1366 endif
1367 amnull=amnull+200. ! ???????
1368 if(ket.ne.0)goto12
1369 endif
1370
1371 c removing bottom mesons to remove b quarks or antiquarks
1372 if(keb.ne.0)then
1373 11 continue
1374 ii=sign(1,keb)
1375 keb=keb-ii
1376 if(ii*keu.le.ii*ked)then
1377 keu=keu+ii
1378 else
1379 ked=ked+ii
1380 endif
1381 amnull=amnull+6. !5.28 ! (more than B-meson)
1382 if(keb.ne.0)goto11
1383 endif
1384
1385 c removing charm mesons to remove c quarks or antiquarks
1386 if(kec.ne.0)then
1387 10 continue
1388 ii=sign(1,kec)
1389 kec=kec-ii
1390 if(keu*ii.le.ked*ii)then
1391 keu=keu+ii
1392 else
1393 ked=ked+ii
1394 endif
1395 amnull=amnull+2.2 !1.87 ! (more than D-meson)
1396 if(kec.ne.0)goto10
1397 endif
1398
1399 c write(ifch,100)keu,ked,kes,kec,keb,ket,amnull
1400
1401 c removing mesons to remove s antiquarks
1402 5 continue
1403 if(kes.lt.0)then
1404 amnull=amnull+asuha(6)
1405 if(keu.ge.ked)then
1406 keu=keu-1
1407 else
1408 ked=ked-1
1409 endif
1410 kes=kes+1
1411 goto5
1412 endif
1413
1414 c removing mesons to remove d antiquarks
1415 6 continue
1416 if(ked.lt.0)then
1417 if(keu.ge.kes)then
1418 amnull=amnull+asuha(5)
1419 keu=keu-1
1420 else
1421 amnull=amnull+asuha(6)
1422 kes=kes-1
1423 endif
1424 ked=ked+1
1425 goto6
1426 endif
1427
1428 c removing mesons to remove u antiquarks
1429 7 continue
1430 if(keu.lt.0)then
1431 if(ked.ge.kes)then
1432 amnull=amnull+asuha(5)
1433 ked=ked-1
1434 else
1435 amnull=amnull+asuha(6)
1436 kes=kes-1
1437 endif
1438 keu=keu+1
1439 goto7
1440 endif
1441
1442 c write(ifch,100)keu,ked,kes,kec,keb,ket,amnull
1443 c print*,keu,ked,kes,kec,keb,ket,amnull
1444
1445 if(keu+ked+kes+kec+keb+ket.ne.ke)
1446 *call utstop('utamnu: sum_kei /= ke&')
1447 keq=keu+ked
1448 keqx=keq
1449 amnux=0
1450
1451 c removing strange baryons
1452 i=4
1453 2 i=i-1
1454 3 continue
1455 if((4-i)*kes.gt.(i-1)*keq)then
1456 amnux=amnux+asuha(1+i)
1457 kes=kes-i
1458 keq=keq-3+i
1459 if(kes.lt.0)call utstop('utamnu: negative kes&')
1460 if(keq.lt.0)call utstop('utamnu: negative keq&')
1461 goto3
1462 endif
1463 if(i.gt.1)goto2
1464 if(keqx.gt.keq)then
1465 do 8 k=1,keqx-keq
1466 if(keu.ge.ked)then
1467 keu=keu-1
1468 else
1469 ked=ked-1
1470 endif
1471 8 continue
1472 endif
1473
1474 if(keu+ked.ne.keq)call utstop('utamnu: keu+ked /= keq&')
1475 c write(ifch,100)keu,ked,kes,kec,keb,ket,amnull+amnux
1476 c print*,keu,ked,kes,kec,keb,ket,amnull+amnux
1477
1478 c removing nonstrange baryons
1479 9 continue
1480 if(keu.gt.2*ked)then
1481 amnux=amnux+asuha(7)
1482 keu=keu-3
1483 if(keu.lt.0)call utstop('utamnu: negative keu&')
1484 goto9
1485 endif
1486 if(ked.gt.2*keu)then
1487 amnux=amnux+asuha(7)
1488 ked=ked-3
1489 if(ked.lt.0)call utstop('utamnu: negative ked&')
1490 goto9
1491 endif
1492 keq=keu+ked
1493
1494 c write(ifch,100)keu,ked,kes,kec,keb,ket,amnull+amnux
1495 c print*,keu,ked,kes,kec,keb,ket,amnull+amnux
1496
1497 if(mod(keq,3).ne.0)call utstop('utamnu: mod(keq,3) /= 0&')
1498 amnux=amnux+asuha(1)*keq/3
1499
1500 c write(ifch,100)keu,ked,kes,kec,keb,ket,amnull+amnux
1501 c print*,keu,ked,kes,kec,keb,ket,amnull+amnux
1502
1503 amnull=amnull+amnux
1504
1505 if(amnull.eq.0)amnull=asuha(5)
1506
1507 utamnu=amnull
1508 return
1509 end
1510
1511 c-----------------------------------------------------------------------
1512 function utamnx(jcp,jcm)
1513 c-----------------------------------------------------------------------
1514 c returns minimum mass for the decay of jcp---jcm (by calling utamnu).
1515 c-----------------------------------------------------------------------
1516 parameter (nflav=6)
1517 integer jcp(nflav,2),jcm(nflav,2)
1518
1519 do i=1,nflav
1520 do j=1,2
1521 if(jcp(i,j).ne.0)goto1
1522 enddo
1523 enddo
1524 keu=jcm(1,1)-jcm(1,2)
1525 ked=jcm(2,1)-jcm(2,2)
1526 kes=jcm(3,1)-jcm(3,2)
1527 kec=jcm(4,1)-jcm(4,2)
1528 keb=jcm(5,1)-jcm(5,2)
1529 ket=jcm(6,1)-jcm(6,2)
1530 utamnx=utamnu(keu,ked,kes,kec,keb,ket,5)
1531 return
1532 1 continue
1533
1534 do i=1,nflav
1535 do j=1,2
1536 if(jcm(i,j).ne.0)goto2
1537 enddo
1538 enddo
1539 keu=jcp(1,1)-jcp(1,2)
1540 ked=jcp(2,1)-jcp(2,2)
1541 kes=jcp(3,1)-jcp(3,2)
1542 kec=jcp(4,1)-jcp(4,2)
1543 keb=jcp(5,1)-jcp(5,2)
1544 ket=jcp(6,1)-jcp(6,2)
1545 utamnx=utamnu(keu,ked,kes,kec,keb,ket,5)
1546 return
1547 2 continue
1548
1549 keu=jcp(1,1)-jcp(1,2)
1550 ked=jcp(2,1)-jcp(2,2)
1551 kes=jcp(3,1)-jcp(3,2)
1552 kec=jcp(4,1)-jcp(4,2)
1553 keb=jcp(5,1)-jcp(5,2)
1554 ket=jcp(6,1)-jcp(6,2)
1555 ke=keu+ked+kes+kec+keb+ket
1556 if(mod(ke+1,3).eq.0)then
1557 keu=keu+1
1558 amms1=utamnu(keu,ked,kes,kec,keb,ket,5)
1559 keu=keu-1
1560 ked=ked+1
1561 amms2=utamnu(keu,ked,kes,kec,keb,ket,5)
1562 elseif(mod(ke-1,3).eq.0)then
1563 keu=keu-1
1564 amms1=utamnu(keu,ked,kes,kec,keb,ket,5)
1565 keu=keu+1
1566 ked=ked-1
1567 amms2=utamnu(keu,ked,kes,kec,keb,ket,5)
1568 else
1569 amms1=0
1570 amms2=0
1571 amms3=0
1572 amms4=0
1573 call utstop('utamnx: no singlet possible (1)&')
1574 endif
1575 keu=jcm(1,1)-jcm(1,2)
1576 ked=jcm(2,1)-jcm(2,2)
1577 kes=jcm(3,1)-jcm(3,2)
1578 kec=jcm(4,1)-jcm(4,2)
1579 keb=jcm(5,1)-jcm(5,2)
1580 ket=jcm(6,1)-jcm(6,2)
1581 ke=keu+ked+kes+kec+keb+ket
1582 if(mod(ke+1,3).eq.0)then
1583 keu=keu+1
1584 amms3=utamnu(keu,ked,kes,kec,keb,ket,5)
1585 keu=keu-1
1586 ked=ked+1
1587 amms4=utamnu(keu,ked,kes,kec,keb,ket,5)
1588 elseif(mod(ke-1,3).eq.0)then
1589 keu=keu-1
1590 amms3=utamnu(keu,ked,kes,kec,keb,ket,5)
1591 keu=keu+1
1592 ked=ked-1
1593 amms4=utamnu(keu,ked,kes,kec,keb,ket,5)
1594 else
1595 call utstop('utamnx: no singlet possible (2)&')
1596 endif
1597 utamnx=min(amms1+amms3,amms2+amms4)
1598 c print *,amms1,amms3,amms2,amms4,jcp,jcm
1599 return
1600 end
1601
1602
1603
1604 cc-----------------------------------------------------------------------
1605 c function utamny(jcp,jcm)
1606 cc-----------------------------------------------------------------------
1607 cc returns minimum mass of jcp+jcm (by calling utamnu).
1608 cc-----------------------------------------------------------------------
1609 c parameter (nflav=6)
1610 c integer jcp(nflav,2),jcm(nflav,2),jc(nflav,2)
1611 c do 7 nf=1,nflav
1612 c jc(nf,1)=jcp(nf,1)+jcm(nf,1)
1613 c7 jc(nf,2)=jcp(nf,2)+jcm(nf,2)
1614 c keu=jc(1,1)-jc(1,2)
1615 c ked=jc(2,1)-jc(2,2)
1616 c kes=jc(3,1)-jc(3,2)
1617 c kec=jc(4,1)-jc(4,2)
1618 c keb=jc(5,1)-jc(5,2)
1619 c ket=jc(6,1)-jc(6,2)
1620 c utamny=utamnu(keu,ked,kes,kec,keb,ket,5)
1621 c return
1622 c end
1623 c
1624 c-----------------------------------------------------------------------
1625 function utamnz(jc,modus)
1626 c-----------------------------------------------------------------------
1627 c returns minimum mass of jc (by calling utamnu).
1628 c-----------------------------------------------------------------------
1629 parameter (nflav=6)
1630 integer jc(nflav,2)
1631 keu=jc(1,1)-jc(1,2)
1632 ked=jc(2,1)-jc(2,2)
1633 kes=jc(3,1)-jc(3,2)
1634 kec=jc(4,1)-jc(4,2)
1635 keb=jc(5,1)-jc(5,2)
1636 ket=jc(6,1)-jc(6,2)
1637 utamnz=utamnu(keu,ked,kes,kec,keb,ket,modus)
1638 return
1639 end
1640
1641 c-----------------------------------------------------------------------
1642 subroutine utar(i1,i2,i3,x0,x1,x2,x3,xx)
1643 c-----------------------------------------------------------------------
1644 c returns the array xx with xx(1)=x0 <= xx(i) <= xx(i3)=x3
1645 c-----------------------------------------------------------------------
1646 real xx(i3)
1647 do 1 i=1,i1-1
1648 1 xx(i)=x0+(i-1.)/(i1-1.)*(x1-x0)
1649 do 2 i=i1,i2-1
1650 2 xx(i)=x1+(i-i1*1.)/(i2-i1*1.)*(x2-x1)
1651 do 3 i=i2,i3
1652 3 xx(i)=x2+(i-i2*1.)/(i3-i2*1.)*(x3-x2)
1653 return
1654 end
1655
1656 cc---------------------------------------------------------------------
1657 c subroutine utaxis(i,j,a1,a2,a3)
1658 cc-----------------------------------------------------------------------
1659 cc calculates the axis defined by the ptls i,j in the i,j cm system
1660 cc---------------------------------------------------------------------
1661 c include 'epos.inc'
1662 c double precision pi1,pi2,pi3,pi4,pj1,pj2,pj3,pj4,p1,p2,p3,p4,p5
1663 c *,err,a
1664 c a1=0
1665 c a2=0
1666 c a3=1
1667 c pi1=dble(pptl(1,i))
1668 c pi2=dble(pptl(2,i))
1669 c pi3=dble(pptl(3,i))
1670 c pi4=dble(pptl(4,i))
1671 c pj1=dble(pptl(1,j))
1672 c pj2=dble(pptl(2,j))
1673 c pj3=dble(pptl(3,j))
1674 c pj4=dble(pptl(4,j))
1675 c p1=pi1+pj1
1676 c p2=pi2+pj2
1677 c p3=pi3+pj3
1678 c p4=pi4+pj4
1679 c p5=dsqrt(p4**2-p3**2-p2**2-p1**2)
1680 c call utlob2(1,p1,p2,p3,p4,p5,pi1,pi2,pi3,pi4,50)
1681 c call utlob2(1,p1,p2,p3,p4,p5,pj1,pj2,pj3,pj4,51)
1682 c err=(pi1+pj1)**2+(pi2+pj2)**2+(pi3+pj3)**2
1683 c if(err.gt.1d-3)then
1684 c call utmsg('utaxis')
1685 c write(ifch,*)'***** err=',err
1686 c write(ifch,*)'pi:',pi1,pi2,pi3,pi4
1687 c write(ifch,*)'pj:',pj1,pj2,pj3,pj4
1688 c call utmsgf
1689 c endif
1690 c a=dsqrt( (pj1-pi1)**2 + (pj2-pi2)**2 + (pj3-pi3)**2 )
1691 c if(a.eq.0.d0)return
1692 c a1=sngl((pi1-pj1)/a)
1693 c a2=sngl((pi2-pj2)/a)
1694 c a3=sngl((pi3-pj3)/a)
1695 c return
1696 c end
1697 c
1698 cc---------------------------------------------------------------------
1699 c subroutine uthalf(i,j,zz1,zz2,iret)
1700 cc-----------------------------------------------------------------------
1701 cc give equal energy (E_i+E_j)/2 to particle i+j in their cm system
1702 cc---------------------------------------------------------------------
1703 c include 'epos.inc'
1704 c double precision pi1,pi2,pi3,pi4,pi5,pj1,pj2,pj3,pj4,pj5
1705 c *,p1,p2,p3,p4,p5,err,pt,pti,sinp,cosp,pmax,phi,drangen!,rrr
1706 c iret=0
1707 c pi1=dble(pptl(1,i))
1708 c pi2=dble(pptl(2,i))
1709 c pi3=dble(pptl(3,i))
1710 c pi5=dble(pptl(5,i))
1711 c pi4=sqrt(pi1**2+pi2**2+pi3**2+pi5**2)
1712 c pj1=dble(pptl(1,j))
1713 c pj2=dble(pptl(2,j))
1714 c pj3=dble(pptl(3,j))
1715 c pj5=dble(pptl(5,j))
1716 c pj4=sqrt(pj1**2+pj2**2+pj3**2+pj5**2)
1717 c if(ish.ge.6)then
1718 c write(ifch,*)'uthalf for ',i,' and ',j
1719 c write(ifch,*)'in ',idptl(i),pi1,pi2,pi3,pi4,pi5
1720 c write(ifch,*)'<> ',idptl(j),pj1,pj2,pj3,pj4,pj5
1721 c endif
1722 c p1=pi1+pj1
1723 c p2=pi2+pj2
1724 c p3=pi3+pj3
1725 c p4=pi4+pj4
1726 c p5=(p4-p3)*(p4+p3)-p2**2-p1**2
1727 c if(p5.lt.0d0.or.(pi3.lt.0.99d0*pj3.and.mod(ityptl(j)/10,10).eq.4)
1728 c & .or.(pi3.gt.0.99d0*pj3.and.mod(ityptl(j)/10,10).eq.5))then
1729 cc if(p5.lt.0d0)then
1730 c if(ish.ge.7)write(ifch,*)'Inversion not possible (1)',p5,pi3,pj3
1731 c iret=1
1732 c return
1733 c else
1734 c p5=sqrt(p5)
1735 c endif
1736 c call utlob2(1,p1,p2,p3,p4,p5,pi1,pi2,pi3,pi4,50)
1737 c call utlob2(1,p1,p2,p3,p4,p5,pj1,pj2,pj3,pj4,51)
1738 c err=(pi1+pj1)**2+(pi2+pj2)**2+(pi3+pj3)**2
1739 c if(err.gt.1d-3)then
1740 c call utmsg('uthalf')
1741 c write(ifch,*)'***** err=',err
1742 c write(ifch,*)'pi:',pi1,pi2,pi3,pi4
1743 c write(ifch,*)'pj:',pj1,pj2,pj3,pj4
1744 c call utmsgf
1745 c endif
1746 c if(ish.ge.8)then
1747 c write(ifch,*)'pi 1:',pi1,pi2,pi3,pi4,pi5
1748 c & ,sqrt(pi1**2+pi2**2+pi3**2+pi5**2)
1749 c write(ifch,*)'pj 1:',pj1,pj2,pj3,pj4,pj5
1750 c & ,sqrt(pj1**2+pj2**2+pj3**2+pj5**2)
1751 c endif
1752 c
1753 c phi=drangen(p5)*2d0*dble(pi)
1754 c pti=sqrt(pi1*pi1+pi2*pi2)
1755 c
1756 cc sinp=abs(asin(pi2/pti))
1757 c
1758 c cosp=cos(phi)
1759 c sinp=sin(phi)
1760 cc cosp=pi1/pti
1761 cc sinp=pi2/pti
1762 c pini=abs(pi3)
1763 c pmax=pini
1764 c
1765 c ntry=0
1766 c 10 ntry=ntry+1
1767 c pi1=-pj1
1768 c pi2=-pj2
1769 c pi3=-pj3
1770 c
1771 c
1772 cc rrr=dble(ranptcut(4.))*max(0.1d0,dble(zz))+0.01d0
1773 cc rrr=dble(ranptcut(zz*max(1.,float(ntry/10)**3))) !flip if pt too large
1774 cc rrr=dble(min(1.,ranptcut(zz/float(max(1,ntry/10))))) !return if pt too large
1775 c
1776 c rrr=rangen()
1777 c phi=dble(pi*zz1*(1.-rrr*zz2/float(max(1,ntry/10))))
1778 c call utroa2(phi,cosp,sinp,0d0,pi1,pi2,pi3) !rotation around an axis perpendicular to p3
1779 c pt=pi1*pi1+pi2*pi2
1780 c
1781 ccc pi3=pini*(1d0-(pmax/pini+1d0)*min(1d0,rrr))
1782 ccc & *sign(1.d0,pi3)
1783 cc pi3=pini*(1d0-(pmax/pini+1d0)*(1d0-min(1d0,rrr)))
1784 cc & *sign(1.d0,pi3)
1785 cccc pi3=(1d0-exp(-drangen(pti)))*pi3
1786 cccc pi3=min(1d0,(drangen()+exp(-min(100d0,dble(zz)))))*pi3
1787 cccc pi3=(dble(reminv)+exp(-min(100d0,dble(zz))))*pi3
1788 cc pt=((pi4+pi3)*(pi4-pi3)-pi5*pi5)
1789 c if(ish.ge.8)write(ifch,*)'ut',ntry,zz,rrr,-pi3/pj3,pt,pti*pti
1790 c &,idptl(i),idptl(j)
1791 c if((pt.lt.0d0.or.pt.gt.max(2.*pti*pti,1.d0))
1792 cc if(pt.lt.0d0
1793 c & .and.ntry.lt.1000)then
1794 c goto 10
1795 c elseif(ntry.ge.1000)then
1796 c if(ish.ge.7)write(ifch,*)'Inversion not possible (2)',pt,pti*pti
1797 c iret=1
1798 c return !pion distribution fall at xf=1
1799 cc pi3=pj3 !if flip all particle with large pt
1800 cc pt=pti*pti !then very very hard pi+ spectra (flat at xf=1 like in pp400 !)
1801 c !but to hard for K, rho, eta, etc ..
1802 c endif
1803 cc print *,'ut',ntry,phi/(pi-sinp),zz,rrr,-pi3/pj3,pt/pti/pti
1804 cc pt=sqrt(pt)
1805 cc pi1=cosp*pt
1806 cc pi2=sinp*pt
1807 c
1808 c pj1=-pi1
1809 c pj2=-pi2
1810 c pj3=-pi3
1811 cc pi4=sqrt(pi3*pi3+pt*pt+pi5*pi5)
1812 cc pj4=sqrt(pj3*pj3+pt*pt+pj5*pj5)
1813 c pi4=sqrt(pi3*pi3+pt+pi5*pi5)
1814 c pj4=sqrt(pj3*pj3+pt+pj5*pj5)
1815 c if(ish.ge.8)then
1816 c write(ifch,*)'pi 2:',pi1,pi2,pi3,pi4,pi5
1817 c & ,sqrt(pi1**2+pi2**2+pi3**2+pi5**2)
1818 c write(ifch,*)'pj 2:',pj1,pj2,pj3,pj4,pj5
1819 c & ,sqrt(pj1**2+pj2**2+pj3**2+pj5**2)
1820 c endif
1821 c call utlob2(-1,p1,p2,p3,p4,p5,pi1,pi2,pi3,pi4,-50)
1822 c call utlob2(-1,p1,p2,p3,p4,p5,pj1,pj2,pj3,pj4,-51)
1823 c if(pi3/dble(pptl(3,i)).gt.1.00001d0
1824 c &.or.(pi3.gt.1.001d0*pj3.and.mod(ityptl(j)/10,10).eq.4)
1825 c &.or.(pi3.lt.1.001d0*pj3.and.mod(ityptl(j)/10,10).eq.5))then
1826 c if(ish.ge.7)write(ifch,*)'Inversion not possible (3)',pi3,pj3
1827 c iret=1
1828 c return
1829 c endif
1830 c id=idptl(i)
1831 c pptl(1,i)=sngl(pi1)
1832 c pptl(2,i)=sngl(pi2)
1833 c pptl(3,i)=sngl(pi3)
1834 c pptl(4,i)=sngl(pi4)
1835 c pptl(5,i)=sngl(pi5)
1836 cc pptl(5,i)=sngl(pj5)
1837 cc idptl(i)=idptl(j)
1838 c pptl(1,j)=sngl(pj1)
1839 c pptl(2,j)=sngl(pj2)
1840 c pptl(3,j)=sngl(pj3)
1841 c pptl(4,j)=sngl(pj4)
1842 c pptl(5,j)=sngl(pj5)
1843 cc pptl(5,j)=sngl(pi5)
1844 cc idptl(j)=id
1845 c if(ish.ge.6)then
1846 c write(ifch,*)'out ',idptl(i),pi1,pi2,pi3,pi4
1847 c write(ifch,*)' <> ',idptl(j),pj1,pj2,pj3,pj4
1848 c endif
1849 c
1850 c return
1851 c end
1852 c
1853 cc-----------------------------------------------------------------------
1854 c subroutine utchm(arp,arm,ii)
1855 cc-----------------------------------------------------------------------
1856 cc checks whether arp**2=0 and arm**2=0.
1857 cc-----------------------------------------------------------------------
1858 c include 'epos.inc'
1859 c double precision arp(4),arm(4),difp,difm
1860 c difp=arp(4)**2-arp(1)**2-arp(2)**2-arp(3)**2
1861 c difm=arm(4)**2-arm(1)**2-arm(2)**2-arm(3)**2
1862 c if(dabs(difp).gt.1e-3*arp(4)**2
1863 c *.or.dabs(difm).gt.1e-3*arm(4)**2)then
1864 c call utmsg('utchm ')
1865 c write(ifch,*)'***** mass non zero - ',ii
1866 c write(ifch,*)'jet-mass**2`s: ',difp,difm
1867 c write(ifch,*)'energy**2`s: ',arp(4)**2,arm(4)**2
1868 c write(ifch,*)(sngl(arp(i)),i=1,4)
1869 c write(ifch,*)(sngl(arm(i)),i=1,4)
1870 c call utmsgf
1871 c endif
1872 c return
1873 c end
1874 c
1875 c-----------------------------------------------------------------------
1876 subroutine utclea(nptlii,nptl0)
1877 c-----------------------------------------------------------------------
1878 c starting from nptlii
1879 c overwrites istptl=99 particles in /cptl/, reduces so nptl
1880 c and update minfra and maxfra
1881 c-----------------------------------------------------------------------
1882 include 'epos.inc'
1883 integer newptl(mxptl)!,oldptl(mxptl),ii(mxptl)
1884
1885 ish0=ish
1886 if(ishsub/100.eq.18)ish=mod(ishsub,100)
1887
1888 call utpri('utclea',ish,ishini,2)
1889
1890 nptli=max(maproj+matarg+1,nptlii)
1891 minfra0=minfra
1892 maxfra0=maxfra
1893 minfra1=maxfra
1894 maxfra1=minfra
1895 if(ish.ge.2)write(ifch,*)'entering subr utclea:',nptl
1896 & ,minfra,maxfra
1897 if(ish.ge.7)then
1898 write(ifch,*)('-',l=1,68)
1899 write(ifch,*)'sr utclea. initial.'
1900 write(ifch,*)('-',l=1,68)
1901 do 34 n=nptli,nptl
1902 write(ifch,116)iorptl(n),jorptl(n),n,ifrptl(1,n),ifrptl(2,n)
1903 *,idptl(n),sqrt(pptl(1,n)**2+pptl(2,n)**2),pptl(3,n),pptl(5,n)
1904 *,istptl(n),ityptl(n)
1905 34 continue
1906 116 format(1x,i6,i6,4x,i6,4x,i6,i6,i12,3x,3(e8.2,1x),i3,i3)
1907 endif
1908
1909 c ish=ish0
1910 c ish0=ish
1911 c if(ishsub/100.eq.18)ish=mod(ishsub,100)
1912
1913 i=nptli-1
1914 1 i=i+1
1915 if(i.gt.nptl)goto 1000
1916 if(istptl(i).eq.99)goto 2
1917 newptl(i)=i
1918 c oldptl(i)=i
1919 goto 1
1920
1921 2 i=i-1
1922 j=i
1923 3 i=i+1
1924 4 j=j+1
1925 if(j.gt.nptl)goto 5
1926 newptl(j)=0
1927 if(istptl(j).eq.99)goto 4
1928 newptl(j)=i
1929 c oldptl(i)=j
1930 c write(ifch,*)'move',j,' to ',i
1931 c write(ifch,*)idptl(i),ityptl(i),idptl(j),ityptl(j),minfra,maxfra
1932 call utrepl(i,j)
1933 if(j.ge.minfra0.and.j.le.maxfra0)then
1934 minfra1=min(minfra1,i)
1935 maxfra1=max(maxfra1,i)
1936 endif
1937 goto 3
1938
1939 5 nptl=i-1
1940 if(nptl.eq.0)then
1941 nptl0=0
1942 goto 1000
1943 endif
1944
1945 20 n0=newptl(nptl0)
1946 if(n0.gt.0)then
1947 nptl0=n0
1948 else
1949 nptl0=nptl0-1
1950 if(nptl0.gt.0)goto 20
1951 endif
1952
1953
1954 c do 11 k=1,nptl
1955 c io=iorptl(k)
1956 c if(io.le.0)ii(k)=io
1957 c if(io.gt.0)ii(k)=newptl(io)
1958 c11 continue
1959 c do 12 k=1,nptl
1960 c12 iorptl(k)=ii(k)
1961 c
1962 c do 13 k=1,nptl
1963 c jo=jorptl(k)
1964 c if(jo.le.0)ii(k)=jo
1965 c if(jo.gt.0)ii(k)=newptl(jo)
1966 c13 continue
1967 c do 14 k=1,nptl
1968 c14 jorptl(k)=ii(k)
1969 c
1970 c do 15 k=1,nptl
1971 c if1=ifrptl(1,k)
1972 c if(if1.le.0)ii(k)=if1
1973 c if(if1.gt.0)ii(k)=newptl(if1)
1974 c15 continue
1975 c do 16 k=1,nptl
1976 c16 ifrptl(1,k)=ii(k)
1977 c
1978 c do 17 k=1,nptl
1979 c if2=ifrptl(2,k)
1980 c if(if2.le.0)ii(k)=if2
1981 c if(if2.gt.0)ii(k)=newptl(if2)
1982 c17 continue
1983 c do 18 k=1,nptl
1984 c18 ifrptl(2,k)=ii(k)
1985 c
1986 c do 19 k=1,nptl
1987 c if(ifrptl(1,k).eq.0.and.ifrptl(2,k).gt.0)ifrptl(1,k)=ifrptl(2,k)
1988 c if(ifrptl(2,k).eq.0.and.ifrptl(1,k).gt.0)ifrptl(2,k)=ifrptl(1,k)
1989 c19 continue
1990
1991 1000 continue
1992
1993 if(minfra1.lt.minfra0)minfra=minfra1
1994 if(maxfra1.ge.minfra1)maxfra=maxfra1
1995
1996 if(ish.ge.2)then
1997 write(ifch,*)'exiting subr utclea:'
1998 do 35 n=1,nptl
1999 write(ifch,116)iorptl(n),jorptl(n),n,ifrptl(1,n),ifrptl(2,n)
2000 *,idptl(n),sqrt(pptl(1,n)**2+pptl(2,n)**2),pptl(3,n),pptl(5,n)
2001 *,istptl(n),ityptl(n)
2002 35 continue
2003 endif
2004
2005 if(ish.ge.2)write(ifch,*)'exiting subr utclea:',nptl
2006 & ,minfra,maxfra
2007
2008 call utprix('utclea',ish,ishini,2)
2009 ish=ish0
2010 return
2011 end
2012
2013 c---------------------------------------------------------------------
2014 subroutine utfit(x,y,ndata,sig,mwt,a,b,siga,sigb,chi2,q)
2015 c---------------------------------------------------------------------
2016 c linear fit to data
2017 c input:
2018 c ndata: nr of data points
2019 c x(),y(),sig(): data
2020 c mwt: unweighted (0) or weighted (else) data points
2021 c output:
2022 c a,b: parameters of linear fit a+b*x
2023 c---------------------------------------------------------------------
2024 INTEGER mwt,ndata
2025 REAL a,b,chi2,q,siga,sigb,sig(ndata),x(ndata),y(ndata)
2026 CU USES utgmq
2027 INTEGER i
2028 REAL sigdat,ss,st2,sx,sxoss,sy,t,wt,utgmq
2029 sx=0.
2030 sy=0.
2031 st2=0.
2032 b=0.
2033 if(mwt.ne.0) then
2034 ss=0.
2035 do 11 i=1,ndata
2036 wt=1./(sig(i)**2)
2037 ss=ss+wt
2038 sx=sx+x(i)*wt
2039 sy=sy+y(i)*wt
2040 11 continue
2041 else
2042 do 12 i=1,ndata
2043 sx=sx+x(i)
2044 sy=sy+y(i)
2045 12 continue
2046 ss=float(ndata)
2047 endif
2048 sxoss=sx/ss
2049 if(mwt.ne.0) then
2050 do 13 i=1,ndata
2051 t=(x(i)-sxoss)/sig(i)
2052 st2=st2+t*t
2053 b=b+t*y(i)/sig(i)
2054 13 continue
2055 else
2056 do 14 i=1,ndata
2057 t=x(i)-sxoss
2058 st2=st2+t*t
2059 b=b+t*y(i)
2060 14 continue
2061 endif
2062 b=b/st2
2063 a=(sy-sx*b)/ss
2064 siga=sqrt((1.+sx*sx/(ss*st2))/ss)
2065 sigb=sqrt(1./st2)
2066 chi2=0.
2067 if(mwt.eq.0) then
2068 do 15 i=1,ndata
2069 chi2=chi2+(y(i)-a-b*x(i))**2
2070 15 continue
2071 q=1.
2072 sigdat=sqrt(chi2/(ndata-2))
2073 siga=siga*sigdat
2074 sigb=sigb*sigdat
2075 else
2076 do 16 i=1,ndata
2077 chi2=chi2+((y(i)-a-b*x(i))/sig(i))**2
2078 16 continue
2079 q=utgmq(0.5*(ndata-2),0.5*chi2)
2080 endif
2081 return
2082 END
2083
2084 c-----------------------------------------------------------------------
2085 function utgam1(x)
2086 c-----------------------------------------------------------------------
2087 c gamma fctn tabulated
2088 c single precision
2089 c-----------------------------------------------------------------------
2090 double precision utgamtab,utgam,al,dl
2091 common/gamtab/utgamtab(10000)
2092
2093 if(x.gt.0.01.and.x.lt.99.99)then
2094 al=100.d0*dble(x)
2095 k1=int(al)
2096 k2=k1+1
2097 dl =al-dble(k1)
2098 utgam1=real(utgamtab(k2)*dl+utgamtab(k1)*(1.d0-dl))
2099 elseif(x.eq.0.)then
2100 utgam1=0.
2101 else
2102 utgam1=real(utgam(dble(x)))
2103 endif
2104
2105 end
2106
2107 c-----------------------------------------------------------------------
2108 double precision function utgam2(x)
2109 c-----------------------------------------------------------------------
2110 c gamma fctn tabulated
2111 c double precision
2112 c-----------------------------------------------------------------------
2113 double precision utgamtab,x,al,dl,utgam
2114 common/gamtab/utgamtab(10000)
2115
2116 if(x.gt.0.01d0.and.x.le.99.99d0)then
2117 al=100.d0*x
2118 k1=int(al)
2119 k2=k1+1
2120 dl =al-dble(k1)
2121 utgam2=utgamtab(k2)*dl+utgamtab(k1)*(1.d0-dl)
2122 elseif(x.eq.0.d0)then
2123 utgam2=0.d0
2124 else
2125 utgam2=utgam(x)
2126 endif
2127
2128 end
2129
2130 c-----------------------------------------------------------------------
2131 double precision function utgam(x)
2132 c-----------------------------------------------------------------------
2133 c gamma fctn
2134 c double precision
2135 c-----------------------------------------------------------------------
2136 include 'epos.inc'
2137 double precision c(13),x,z,f
2138 data c
2139 1/ 0.00053 96989 58808, 0.00261 93072 82746, 0.02044 96308 23590,
2140 2 0.07309 48364 14370, 0.27964 36915 78538, 0.55338 76923 85769,
2141 3 0.99999 99999 99998,-0.00083 27247 08684, 0.00469 86580 79622,
2142 4 0.02252 38347 47260,-0.17044 79328 74746,-0.05681 03350 86194,
2143 5 1.13060 33572 86556/
2144 utgam=0d0
2145 z=x
2146 if(x .gt. 170.d0) goto6
2147 if(x .gt. 0.0d0) goto1
2148 if(x .eq. int(x)) goto5
2149 z=1.0d0-z
2150 1 f=1.0d0/z
2151 if(z .le. 1.0d0) goto4
2152 f=1.0d0
2153 2 continue
2154 if(z .lt. 2.0d0) goto3
2155 z=z-1.0d0
2156 f=f*z
2157 goto2
2158 3 z=z-1.0d0
2159 4 utgam=
2160 1 f*((((((c(1)*z+c(2))*z+c(3))*z+c(4))*z+c(5))*z+c(6))*z+c(7))/
2161 2 ((((((c(8)*z+c(9))*z+c(10))*z+c(11))*z+c(12))*z+c(13))*z+1.0d0)
2162 if(x .gt. 0.0d0) return
2163 utgam=3.141592653589793d0/(sin(3.141592653589793d0*x)*utgam)
2164 return
2165 5 write(ifch,10)sngl(x)
2166 10 format(1x,'argument of gamma fctn = ',e20.5)
2167 call utstop('utgam : negative integer argument&')
2168 6 write(ifch,11)sngl(x)
2169 11 format(1x,'argument of gamma fctn = ',e20.5)
2170 call utstop('utgam : argument too large&')
2171 end
2172
2173 c---------------------------------------------------------------------
2174 subroutine utgcf(gammcf,a,x,gln)
2175 c---------------------------------------------------------------------
2176 INTEGER ITMAX
2177 REAL a,gammcf,gln,x,EPS,FPMIN
2178 PARAMETER (ITMAX=100,EPS=3.e-7,FPMIN=1.e-30)
2179 CU USES utgmln
2180 INTEGER i
2181 REAL an,b,c,d,del,h,utgmln
2182 gln=utgmln(a)
2183 b=x+1.-a
2184 c=1./FPMIN
2185 d=1./b
2186 h=d
2187 do 11 i=1,ITMAX
2188 an=-i*(i-a)
2189 b=b+2.
2190 d=an*d+b
2191 if(abs(d).lt.FPMIN)d=FPMIN
2192 c=b+an/c
2193 if(abs(c).lt.FPMIN)c=FPMIN
2194 d=1./d
2195 del=d*c
2196 h=h*del
2197 if(abs(del-1.).lt.EPS)goto 1
2198 11 continue
2199 call utstop("a too large, ITMAX too small in utgcf&")
2200 1 gammcf=exp(-x+a*log(x)-gln)*h
2201 return
2202 END
2203
2204 c---------------------------------------------------------------------
2205 function utgmln(xx)
2206 c---------------------------------------------------------------------
2207 REAL utgmln,xx
2208 INTEGER j
2209 DOUBLE PRECISION ser,stp,tmp,x,y,cof(6)
2210 SAVE cof,stp
2211 DATA cof,stp/76.18009172947146d0,-86.50532032941677d0,
2212 *24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,
2213 *-.5395239384953d-5,2.5066282746310005d0/
2214 x=xx
2215 y=x
2216 tmp=x+5.5d0
2217 tmp=(x+0.5d0)*log(tmp)-tmp
2218 ser=1.000000000190015d0
2219 do 11 j=1,6
2220 y=y+1.d0
2221 ser=ser+cof(j)/y
2222 11 continue
2223 utgmln=tmp+log(stp*ser/x)
2224 return
2225 END
2226
2227 c---------------------------------------------------------------------
2228 function utgmq(a,x)
2229 c---------------------------------------------------------------------
2230 REAL a,utgmq,x
2231 CU USES utgcf,utgser
2232 REAL gammcf,gamser,gln
2233 if(x.lt.0..or.a.le.0.) call utstop("bad arguments in utgmq&")
2234 if(x.lt.a+1.)then
2235 call utgser(gamser,a,x,gln)
2236 utgmq=1.-gamser
2237 else
2238 call utgcf(gammcf,a,x,gln)
2239 utgmq=gammcf
2240 endif
2241 return
2242 END
2243
2244 c---------------------------------------------------------------------
2245 subroutine utgser(gamser,a,x,gln)
2246 c---------------------------------------------------------------------
2247 INTEGER ITMAX
2248 REAL a,gamser,gln,x,EPS
2249 PARAMETER (ITMAX=100,EPS=3.e-7)
2250 CU USES utgmln
2251 INTEGER n
2252 REAL ap,del,sum,utgmln
2253 gln=utgmln(a)
2254 if(x.le.0.)then
2255 if(x.lt.0.)call utstop("x < 0 in utgser&")
2256 gamser=0.
2257 return
2258 endif
2259 ap=a
2260 sum=1./a
2261 del=sum
2262 do 11 n=1,ITMAX
2263 ap=ap+1.
2264 del=del*x/ap
2265 sum=sum+del
2266 if(abs(del).lt.abs(sum)*EPS)goto 1
2267 11 continue
2268 call utstop("a too large, ITMAX too small in utgser&")
2269 1 gamser=sum*exp(-x+a*log(x)-gln)
2270 return
2271 END
2272
2273 c-------------------------------------------------------------------------
2274 subroutine uticpl(ic,ifla,iqaq,iret)
2275 c-------------------------------------------------------------------------
2276 c adds a quark (iqaq=1) or antiquark (iqaq=2) of flavour ifla
2277 c to 2-id ic
2278 c-------------------------------------------------------------------------
2279 include 'epos.inc'
2280 integer jc(nflav,2),ic(2)
2281 iret=0
2282 if(ifla.eq.0)return
2283 call iddeco(ic,jc)
2284 if(ish.ge.8)write(ifch,'(2i8,12i3)')ic,jc
2285 jqaq=3-iqaq
2286 if(jc(ifla,jqaq).gt.0)then
2287 jc(ifla,jqaq)=jc(ifla,jqaq)-1
2288 else
2289 jc(ifla,iqaq)=jc(ifla,iqaq)+1
2290 endif
2291 call idcomj(jc)
2292 call idenco(jc,ic,ireten)
2293 if(ish.ge.8)write(ifch,'(2i8,12i3)')ic,jc
2294 if(ireten.eq.1)iret=1
2295 if(ic(1).eq.0.and.ic(2).eq.0.and.ireten.eq.0)then
2296 ic(1)=100000
2297 ic(2)=100000
2298 endif
2299 return
2300 end
2301
2302 cc-----------------------------------------------------------------------
2303 c subroutine utindx(n,xar,x,i)
2304 cc-----------------------------------------------------------------------
2305 cc input: dimension n
2306 cc array xar(n) with xar(i) > xar(i-1)
2307 cc some number x between xar(1) and xar(n)
2308 cc output: the index i such that x is between xar(i) and xar(i+1)
2309 cc-----------------------------------------------------------------------
2310 c include 'epos.inc'
2311 c real xar(n)
2312 c if(x.lt.xar(1))then
2313 c if(ish.ge.5)then
2314 c call utmsg('utindx')
2315 c write(ifch,*)'***** x=',x,' < xar(1)=',xar(1)
2316 c call utmsgf
2317 c endif
2318 c i=1
2319 c return
2320 c elseif(x.gt.xar(n))then
2321 c if(ish.ge.5)then
2322 c call utmsg('utindx')
2323 c write(ifch,*)'***** x=',x,' > xar(n)=',xar(n)
2324 c call utmsgf
2325 c endif
2326 c i=n
2327 c return
2328 c endif
2329 c lu=1
2330 c lo=n
2331 c1 lz=(lo+lu)/2
2332 c if((xar(lu).le.x).and.(x.le.xar(lz)))then
2333 c lo=lz
2334 c elseif((xar(lz).lt.x).and.(x.le.xar(lo)))then
2335 c lu=lz
2336 c else
2337 c call utstop('utindx: no interval found&')
2338 c endif
2339 c if((lo-lu).ge.2) goto1
2340 c if(lo.le.lu)call utstop('utinvt: lo.le.lu&')
2341 c i=lu
2342 c return
2343 c end
2344 c
2345 c-----------------------------------------------------------------------
2346 function utinvt(n,x,q,y)
2347 c-----------------------------------------------------------------------
2348 c returns x with y=q(x)
2349 c-----------------------------------------------------------------------
2350 include 'epos.inc'
2351 real x(n),q(n)
2352 if(q(n).eq.0.)call utstop('utinvt: q(n)=0&')
2353 if(y.lt.0.)then
2354 if(ish.ge.1)then
2355 call utmsg('utinvt')
2356 write(ifch,*)'***** y=',y,' < 0'
2357 call utmsgf
2358 endif
2359 y=0.
2360 elseif(y.gt.q(n))then
2361 if(ish.ge.1)then
2362 call utmsg('utinvt')
2363 write(ifch,*)'***** y=',y,' > ',q(n)
2364 call utmsgf
2365 endif
2366 y=q(n)
2367 endif
2368 lu=1
2369 lo=n
2370 1 lz=(lo+lu)/2
2371 if((q(lu).le.y).and.(y.le.q(lz)))then
2372 lo=lz
2373 elseif((q(lz).lt.y).and.(y.le.q(lo)))then
2374 lu=lz
2375 else
2376 write(ifch,*)'q(1),y,q(n):',q(1),y,q(n)
2377 write(ifch,*)'lu,lz,lo:',lu,lz,lo
2378 write(ifch,*)'q(lu),q(lz),q(lo):',q(lu),q(lz),q(lo)
2379 call utstop('utinvt: no interval found&')
2380 endif
2381 if((lo-lu).ge.2) goto1
2382 if(lo.le.lu)call utstop('utinvt: lo.le.lu&')
2383 utinvt=x(lu)+(y-q(lu))*(x(lo)-x(lu))/(q(lo)-q(lu))
2384 return
2385 end
2386
2387 c-----------------------------------------------------------------------
2388 subroutine utlob2(isig,p1,p2,p3,p4,p5,x1,x2,x3,x4,idi)
2389 c-----------------------------------------------------------------------
2390 c performs a lorentz boost, double prec.
2391 c isig=+1 is to boost the four vector x1,x2,x3,x4 such as to obtain it
2392 c in the frame specified by the 5-vector p1...p5 (5-vector=4-vector+mass).
2393 c isig=-1: the other way round, that means,
2394 c if the 4-vector x1...x4 is given in some frame characterized by
2395 c p1...p5 with respect to to some lab-frame, utlob2 returns the 4-vector
2396 c x1...x4 in the lab frame.
2397 c idi is a call identifyer (integer) to identify the call in case of problem
2398 c-----------------------------------------------------------------------
2399 include 'epos.inc'
2400 double precision beta(4),z(4),p1,p2,p3,p4,p5,pp,bp,x1,x2,x3,x4
2401 *,xx0,x10,x20,x30,x40,x4x,x0123
2402 if(ish.ge.2)then
2403 if(ish.ge.9)then
2404 write(ifch,101)x1,x2,x3,x4,(x4-x3)*(x4+x3)-x2*x2-x1*x1
2405 write(ifch,301)p1,p2,p3,p4,p5,(p4-p3)*(p4+p3)-p2*p2-p1*p1
2406 101 format(' utlob2: x = ',5e13.5)
2407 301 format(' p = ',6e13.5)
2408 endif
2409 pp=(p4-p3)*(p4+p3)-p2*p2-p1*p1
2410 if(dabs(pp-p5*p5).gt.1e-3*p4*p4.and.dabs(pp-p5*p5).gt.1e-3)then
2411 call utmsg('utlob2')
2412 write(ifch,*)'***** p**2 .ne. p5**2'
2413 write(ifch,*)'call identifyer:',idi
2414 write(ifch,*)'p**2,p5**2: ',pp,p5*p5
2415 write(ifch,*)'p: ',p1,p2,p3,p4,p5
2416 call utmsgf
2417 endif
2418 x10=x1
2419 x20=x2
2420 x30=x3
2421 x40=x4
2422 endif
2423 xx0=(x4-x3)*(x4+x3)-x2*x2-x1*x1
2424 if(p5.le.0.)then
2425 call utmsg('utlob2')
2426 write(ifch,*)'***** p5 negative.'
2427 write(ifch,*)'call identifyer:',idi
2428 write(ifch,*)'p(5): ',p1,p2,p3,p4,p5
2429 write(ifmt,*)'call identifyer:',idi
2430 write(ifmt,*)'p(5): ',p1,p2,p3,p4,p5
2431 call utmsgf
2432 call utstop('utlob2: p5 negative.&')
2433 endif
2434 z(1)=x1
2435 z(2)=x2
2436 z(3)=x3
2437 z(4)=x4
2438 beta(1)=-p1/p5
2439 beta(2)=-p2/p5
2440 beta(3)=-p3/p5
2441 beta(4)= p4/p5
2442 bp=0.
2443 do 220 k=1,3
2444 220 bp=bp+z(k)*isig*beta(k)
2445 do 230 k=1,3
2446 230 z(k)=z(k)+isig*beta(k)*z(4)
2447 *+isig*beta(k)*bp/(beta(4)+1.)
2448 z(4)=beta(4)*z(4)+bp
2449 x1=z(1)
2450 x2=z(2)
2451 x3=z(3)
2452 x4=z(4)
2453 if(ish.ge.9)
2454 *write(ifch,101)x1,x2,x3,x4,(x4-x3)*(x4+x3)-x2*x2-x1*x1
2455 x4x=x4
2456 x0123=xx0+x1*x1+x2*x2+x3*x3
2457 if(x0123.gt.0.)then
2458 x4=sign( dsqrt(x0123) , x4x )
2459 else
2460 x4=0
2461 endif
2462 if(ish.ge.9)then
2463 write(ifch,101)x1,x2,x3,x4,(x4-x3)*(x4+x3)-x2*x2-x1*x1
2464 endif
2465 if(ish.ge.2)then
2466 if(ish.ge.9)write(ifch,*)'check x**2_ini -- x**2_fin'
2467 if(dabs(x4-x4x).gt.1d-2*dabs(x4).and.dabs(x4-x4x).gt.1d-2)then
2468 call utmsg('utlob2')
2469 write(ifch,*)'***** x**2_ini .ne. x**2_fin.'
2470 write(ifch,*)'call identifyer:',idi
2471 write(ifch,*)'x1 x2 x3 x4 x**2 (initial/final/corrected):'
2472 102 format(5e13.5)
2473 write(ifch,102)x10,x20,x30,x40,(x40-x30)*(x40+x30)-x20*x20-x10*x10
2474 write(ifch,102)x1,x2,x3,x4x,(x4x-x3)*(x4x+x3)-x2*x2-x1*x1
2475 write(ifch,102)x1,x2,x3,x4,(x4-x3)*(x4+x3)-x2*x2-x1*x1
2476 call utmsgf
2477 endif
2478 endif
2479 if(ish.ge.9)write(ifch,*)'return from utlob2'
2480 return
2481 end
2482
2483 c-----------------------------------------------------------------------
2484 subroutine utlob3(isig,p1,p2,p3,p4,p5,x1,x2,x3,x4)
2485 c-----------------------------------------------------------------------
2486 c performs a lorentz boost, double prec.
2487 c but arguments are single precision
2488 c-----------------------------------------------------------------------
2489 double precision xx1,xx2,xx3,xx4
2490 xx1=dble(x1)
2491 xx2=dble(x2)
2492 xx3=dble(x3)
2493 xx4=dble(x4)
2494 call utlob2(isig
2495 *,dble(p1),dble(p2),dble(p3),dble(p4),dble(p5)
2496 *,xx1,xx2,xx3,xx4,52)
2497 x1=sngl(xx1)
2498 x2=sngl(xx2)
2499 x3=sngl(xx3)
2500 x4=sngl(xx4)
2501 return
2502 end
2503
2504 c-----------------------------------------------------------------------
2505 subroutine utlob5(yboost,x1,x2,x3,x4,x5)
2506 c-----------------------------------------------------------------------
2507 amt=sqrt(x5**2+x1**2+x2**2)
2508 y=sign(1.,x3)*alog((x4+abs(x3))/amt)
2509 y=y-yboost
2510 x4=amt*cosh(y)
2511 x3=amt*sinh(y)
2512 return
2513 end
2514
2515 c-----------------------------------------------------------------------
2516 subroutine utlob4(isig,pp1,pp2,pp3,pp4,pp5,x1,x2,x3,x4)
2517 c-----------------------------------------------------------------------
2518 c performs a lorentz boost, double prec.
2519 c but arguments are partly single precision
2520 c-----------------------------------------------------------------------
2521 double precision xx1,xx2,xx3,xx4,pp1,pp2,pp3,pp4,pp5
2522 xx1=dble(x1)
2523 xx2=dble(x2)
2524 xx3=dble(x3)
2525 xx4=dble(x4)
2526 call utlob2(isig,pp1,pp2,pp3,pp4,pp5,xx1,xx2,xx3,xx4,53)
2527 x1=sngl(xx1)
2528 x2=sngl(xx2)
2529 x3=sngl(xx3)
2530 x4=sngl(xx4)
2531 return
2532 end
2533
2534
2535 c-----------------------------------------------------------------------
2536 subroutine utlobo(isig,p1,p2,p3,p4,p5,x1,x2,x3,x4)
2537 c-----------------------------------------------------------------------
2538 c performs a lorentz boost
2539 c-----------------------------------------------------------------------
2540 include 'epos.inc'
2541 real beta(4),z(4)
2542 if(p5.le.0.)then
2543 call utmsg('utlobo')
2544 write(ifch,*)'***** mass <= 0.'
2545 write(ifch,*)'p(5): ',p1,p2,p3,p4,p5
2546 call utmsgf
2547 call utstop('utlobo: mass <= 0.&')
2548 endif
2549 z(1)=x1
2550 z(2)=x2
2551 z(3)=x3
2552 z(4)=x4
2553 beta(1)=-p1/p5
2554 beta(2)=-p2/p5
2555 beta(3)=-p3/p5
2556 beta(4)= p4/p5
2557 bp=0.
2558 do 220 k=1,3
2559 220 bp=bp+z(k)*isig*beta(k)
2560 do 230 k=1,3
2561 230 z(k)=z(k)+isig*beta(k)*z(4)
2562 *+isig*beta(k)*bp/(beta(4)+1.)
2563 z(4)=beta(4)*z(4)+bp
2564 x1=z(1)
2565 x2=z(2)
2566 x3=z(3)
2567 x4=z(4)
2568 return
2569 end
2570
2571 c-----------------------------------------------------------------------
2572 subroutine utloc(ar,n,a,l)
2573 c-----------------------------------------------------------------------
2574 real ar(n)
2575 do 1 i=1,n
2576 l=i-1
2577 if(a.lt.ar(i))return
2578 1 continue
2579 l=n
2580 return
2581 end
2582
2583 cc-----------------------------------------------------------------------
2584 c subroutine utlow(cone)
2585 cc-----------------------------------------------------------------------
2586 c character*1 cone
2587 c if(cone.eq.'A')cone='a'
2588 c if(cone.eq.'B')cone='b'
2589 c if(cone.eq.'C')cone='c'
2590 c if(cone.eq.'D')cone='d'
2591 c if(cone.eq.'E')cone='e'
2592 c if(cone.eq.'F')cone='f'
2593 c if(cone.eq.'G')cone='g'
2594 c if(cone.eq.'H')cone='h'
2595 c if(cone.eq.'I')cone='i'
2596 c if(cone.eq.'J')cone='j'
2597 c if(cone.eq.'K')cone='k'
2598 c if(cone.eq.'L')cone='l'
2599 c if(cone.eq.'M')cone='m'
2600 c if(cone.eq.'N')cone='n'
2601 c if(cone.eq.'O')cone='o'
2602 c if(cone.eq.'P')cone='p'
2603 c if(cone.eq.'Q')cone='q'
2604 c if(cone.eq.'R')cone='r'
2605 c if(cone.eq.'S')cone='s'
2606 c if(cone.eq.'T')cone='t'
2607 c if(cone.eq.'U')cone='u'
2608 c if(cone.eq.'V')cone='v'
2609 c if(cone.eq.'W')cone='w'
2610 c if(cone.eq.'X')cone='x'
2611 c if(cone.eq.'Y')cone='y'
2612 c if(cone.eq.'Z')cone='z'
2613 c return
2614 c end
2615 c
2616 cc-----------------------------------------------------------------------
2617 c subroutine utlow3(cthree)
2618 cc-----------------------------------------------------------------------
2619 c character cthree*3
2620 c do 1 i=1,3
2621 c1 call utlow(cthree(i:i))
2622 c return
2623 c end
2624 c
2625 cc-----------------------------------------------------------------------
2626 c subroutine utlow6(csix)
2627 cc-----------------------------------------------------------------------
2628 c character csix*6
2629 c do 1 i=1,6
2630 c1 call utlow(csix(i:i))
2631 c return
2632 c end
2633 c
2634 cc-----------------------------------------------------------------------
2635 c function utmom(k,n,x,q)
2636 cc-----------------------------------------------------------------------
2637 cc calculates kth moment for f(x) with q(i)=int[0,x(i)]f(z)dz
2638 cc-----------------------------------------------------------------------
2639 c real x(n),q(n)
2640 c if(n.lt.2)call utstop('utmom : dimension too small&')
2641 c utmom=0
2642 c do 1 i=2,n
2643 c1 utmom=utmom+((x(i)+x(i-1))/2)**k*(q(i)-q(i-1))
2644 c utmom=utmom/q(n)
2645 c return
2646 c end
2647 c
2648 c-----------------------------------------------------------------------
2649 function utpcm(a,b,c)
2650 c-----------------------------------------------------------------------
2651 c calculates cm momentum for a-->b+c
2652 c-----------------------------------------------------------------------
2653 val=(a*a-b*b-c*c)*(a*a-b*b-c*c)-(2.*b*c)*(2.*b*c)
2654 if(val.lt.0..and.val.gt.-1e-4)then
2655 utpcm=0
2656 return
2657 endif
2658 utpcm=sqrt(val)/(2.*a)
2659 return
2660 end
2661
2662 c-----------------------------------------------------------------------
2663 double precision function utpcmd(a,b,c,iret)
2664 c-----------------------------------------------------------------------
2665 c calculates cm momentum for a-->b+c
2666 c-----------------------------------------------------------------------
2667 double precision a,b,c,val
2668 iret=0
2669 val=(a*a-b*b-c*c)*(a*a-b*b-c*c)-(2.*b*c)*(2.*b*c)
2670 utpcmd=0d0
2671 if(val.lt.0d0.and.val.gt.-1d-4)then
2672 return
2673 elseif(val.lt.0d0)then
2674 iret=1
2675 return
2676 endif
2677 utpcmd=sqrt(val)/(2.d0*a)
2678 return
2679 end
2680
2681 c-----------------------------------------------------------------------
2682 subroutine utpri(text,ishi,ishini,ishx)
2683 c-----------------------------------------------------------------------
2684 include 'epos.inc'
2685 character*6 text
2686 c double precision seedx !!!
2687 ishini=ishi
2688 if(ishevt.ne.0.and.nrevt+1.ne.ishevt)return
2689 if(nrpri.gt.0)then
2690 do nr=1,nrpri
2691 if(subpri(nr)(1:6).eq.text)then
2692 ishi=ishpri(nr)
2693 endif
2694 enddo
2695 endif
2696 if(ish.ge.ishx)then
2697 write(ifch,'(1x,43a)')
2698 * ('-',i=1,10),' entry ',text,' ',('-',i=1,30)
2699 c call ranfgt(seedx) !!!
2700 c if(ish.ge.ishx)write(ifch,*)'seed:',seedx !!!
2701 endif
2702 return
2703 end
2704
2705 c-----------------------------------------------------------------------
2706 subroutine utprix(text,ishi,ishini,ishx)
2707 c-----------------------------------------------------------------------
2708 include 'epos.inc'
2709 character*6 text
2710 if(ishevt.ne.0.and.nrevt+1.ne.ishevt)return
2711 if(ish.ge.ishx)write(ifch,'(1x,44a)')
2712 *('-',i=1,30),' exit ',text,' ',('-',i=1,11)
2713 ishi=ishini
2714 return
2715 end
2716
2717 c-----------------------------------------------------------------------
2718 subroutine utprj(text,ishi,ishini,ishx)
2719 c-----------------------------------------------------------------------
2720 include 'epos.inc'
2721 character*20 text
2722 c double precision seedx !!!
2723 idx=index(text,' ')-1
2724 ishini=ishi
2725 if(ishevt.ne.0.and.nrevt+1.ne.ishevt)return
2726 if(nrpri.gt.0)then
2727 do nr=1,nrpri
2728 if(subpri(nr)(1:idx).eq.text(1:idx))then
2729 ishi=ishpri(nr)
2730 endif
2731 enddo
2732 endif
2733 if(ish.ge.ishx)then
2734 write(ifch,'(1x,43a)')
2735 * ('-',i=1,10),' entry ',text(1:idx),' ',('-',i=1,30)
2736 c call ranfgt(seedx) !!!
2737 c if(ish.ge.ishx)write(ifch,*)'seed:',seedx !!!
2738 endif
2739 return
2740 end
2741
2742 c-----------------------------------------------------------------------
2743 subroutine utprjx(text,ishi,ishini,ishx)
2744 c-----------------------------------------------------------------------
2745 include 'epos.inc'
2746 character*20 text
2747 idx=index(text,' ')-1
2748 if(ishevt.ne.0.and.nrevt+1.ne.ishevt)return
2749 if(ish.ge.ishx)write(ifch,'(1x,44a)')
2750 *('-',i=1,30),' exit ',text(1:idx),' ',('-',i=1,11)
2751 ishi=ishini
2752 return
2753 end
2754
2755 c-----------------------------------------------------------------------
2756 function utquad(m,x,f,k)
2757 c-----------------------------------------------------------------------
2758 c performs an integration according to simpson
2759 c-----------------------------------------------------------------------
2760 real x(m),f(m)
2761 utquad=0
2762 do 1 i=1,k-1
2763 1 utquad=utquad+(f(i)+f(i+1))/2*(x(i+1)-x(i))
2764 return
2765 end
2766
2767 c-----------------------------------------------------------------------
2768 subroutine utquaf(fu,n,x,q,x0,x1,x2,x3)
2769 c-----------------------------------------------------------------------
2770 c returns q(i) = integral [x(1)->x(i)] fu(x) dx
2771 c-----------------------------------------------------------------------
2772 include 'epos.inc'
2773 real x(n),q(n)
2774 parameter (m=10)
2775 real xa(m),fa(m)
2776 external fu
2777 if(x1.lt.x0.or.x2.lt.x1.or.x3.lt.x2)then
2778 if(ish.ge.1)then
2779 call utmsg('utquaf')
2780 write(ifch,*)' xi=',x0,x1,x2,x3
2781 call utmsgf
2782 endif
2783 endif
2784 call utar(n/3,n*2/3,n,x0,x1,x2,x3,x)
2785 q(1)=0
2786 do 2 i=2,n
2787 do 3 k=1,m
2788 z=x(i-1)+(k-1.)/(m-1.)*(x(i)-x(i-1))
2789 xa(k)=z
2790 3 fa(k)=fu(z)
2791 q(i)=q(i-1)+utquad(m,xa,fa,m)
2792 2 continue
2793 return
2794 end
2795
2796 c-----------------------------------------------------------------------
2797 subroutine utrepl(i,j)
2798 c-----------------------------------------------------------------------
2799 c i is replaced by j in /cptl/
2800 c-----------------------------------------------------------------------
2801 include 'epos.inc'
2802 do 1 k=1,5
2803 1 pptl(k,i) =pptl(k,j)
2804 iorptl(i) = 0 !iorptl(j)
2805 idptl(i) =idptl(j)
2806 istptl(i) =istptl(j)
2807 do 2 k=1,2
2808 2 tivptl(k,i)=tivptl(k,j)
2809 do 3 k=1,2
2810 3 ifrptl(k,i)= 0 !ifrptl(k,j)
2811 jorptl(i) = 0 !jorptl(j)
2812 do 4 k=1,4
2813 4 xorptl(k,i)=xorptl(k,j)
2814 do 5 k=1,4
2815 5 ibptl(k,i) =ibptl(k,j)
2816 ityptl(i) =ityptl(j)
2817 iaaptl(i) =iaaptl(j)
2818 radptl(i) =radptl(j)
2819 desptl(i) =desptl(j)
2820 dezptl(i) =dezptl(j)
2821 qsqptl(i) =qsqptl(j)
2822 zpaptl(1,i)=zpaptl(1,j)
2823 zpaptl(2,i)=zpaptl(2,j)
2824 itsptl(i) =itsptl(j)
2825 rinptl(i) =rinptl(j)
2826 return
2827 end
2828
2829 c-----------------------------------------------------------------------
2830 subroutine utrepla(i,j)
2831 c-----------------------------------------------------------------------
2832 c i is replaced by j in /cptl/
2833 c-----------------------------------------------------------------------
2834 include 'epos.inc'
2835 do 1 k=1,5
2836 1 pptl(k,i) =pptl(k,j)
2837 iorptl(i) = iorptl(j)
2838 idptl(i) =idptl(j)
2839 istptl(i) =istptl(j)
2840 do 2 k=1,2
2841 2 tivptl(k,i)=tivptl(k,j)
2842 do 3 k=1,2
2843 3 ifrptl(k,i)= ifrptl(k,j)
2844 jorptl(i) = jorptl(j)
2845 do 4 k=1,4
2846 4 xorptl(k,i)=xorptl(k,j)
2847 do 5 k=1,4
2848 5 ibptl(k,i) =ibptl(k,j)
2849 ityptl(i) =ityptl(j)
2850 iaaptl(i) =iaaptl(j)
2851 radptl(i) =radptl(j)
2852 desptl(i) =desptl(j)
2853 dezptl(i) =dezptl(j)
2854 qsqptl(i) =qsqptl(j)
2855 zpaptl(1,i)=zpaptl(1,j)
2856 zpaptl(2,i)=zpaptl(2,j)
2857 itsptl(i) =itsptl(j)
2858 rinptl(i) =rinptl(j)
2859 return
2860 end
2861
2862 cc-----------------------------------------------------------------------
2863 c subroutine utresm(icp1,icp2,icm1,icm2,amp,idpr,iadj,ireten)
2864 cc-----------------------------------------------------------------------
2865 c parameter (nflav=6)
2866 c integer icm(2),icp(2),jcm(nflav,2),jcp(nflav,2)
2867 c icm(1)=icm1
2868 c icm(2)=icm2
2869 c icp(1)=icp1
2870 c icp(2)=icp2
2871 c CALL IDDECO(ICM,JCM)
2872 c CALL IDDECO(ICP,JCP)
2873 c do 37 nf=1,nflav
2874 c do 37 k=1,2
2875 c37 jcP(nf,k)=jcp(nf,k)+jcm(nf,k)
2876 c CALL IDENCO(JCP,ICP,IRETEN)
2877 c IDP=IDTRA(ICP,0,0,3)
2878 c call idres(idp,amp,idpr,iadj)
2879 c return
2880 c end
2881 c
2882 cc-----------------------------------------------------------------------
2883 c subroutine utroa1(phi,a1,a2,a3,x1,x2,x3)
2884 cc-----------------------------------------------------------------------
2885 cc rotates x by angle phi around axis a.
2886 cc normalization of a is irrelevant.
2887 cc-----------------------------------------------------------------------
2888 c double precision aaa,aa(3),xxx,xx(3),e1(3),e2(3),e3(3),xp,xt,dphi
2889 c dphi=phi
2890 c xx(1)=x1
2891 c xx(2)=x2
2892 c xx(3)=x3
2893 c aa(1)=a1
2894 c aa(2)=a2
2895 c aa(3)=a3
2896 c aaa=0
2897 c xxx=0
2898 c do i=1,3
2899 c aaa=aaa+aa(i)**2
2900 c xxx=xxx+xx(i)**2
2901 c enddo
2902 c if(xxx.eq.0d0)return
2903 c if(aaa.eq.0d0)call utstop('utroa1: zero rotation axis&')
2904 c aaa=dsqrt(aaa)
2905 c xxx=dsqrt(xxx)
2906 cc e3 = a / !a!
2907 c do i=1,3
2908 c e3(i)=aa(i)/aaa
2909 c enddo
2910 cc x_parallel
2911 c xp=0
2912 c do i=1,3
2913 c xp=xp+xx(i)*e3(i)
2914 c enddo
2915 cc x_transverse
2916 c if(xxx**2-xp**2.le.0.)return
2917 c xt=dsqrt(xxx**2-xp**2)
2918 cc e1 = vector x_transverse / absolute value x_transverse
2919 c do i=1,3
2920 c e1(i)=(xx(i)-e3(i)*xp)/xt
2921 c enddo
2922 cc e2 orthogonal e3,e1
2923 c call utvec2(e3,e1,e2)
2924 cc rotate x
2925 c do i=1,3
2926 c xx(i)=xp*e3(i)+xt*dcos(dphi)*e1(i)+xt*dsin(dphi)*e2(i)
2927 c enddo
2928 cc back to single precision
2929 c x1=xx(1)
2930 c x2=xx(2)
2931 c x3=xx(3)
2932 c return
2933 c end
2934 c
2935 c-----------------------------------------------------------------------
2936 subroutine utroa1(phi,a1,a2,a3,x1,x2,x3)
2937 c-----------------------------------------------------------------------
2938 c rotates x by angle phi around axis a (argument single precision)
2939 c normalization of a is irrelevant.
2940 c-----------------------------------------------------------------------
2941 double precision aa(3),xx(3),dphi
2942 dphi=phi
2943 xx(1)=x1
2944 xx(2)=x2
2945 xx(3)=x3
2946 aa(1)=a1
2947 aa(2)=a2
2948 aa(3)=a3
2949 call utroa2(dphi,aa(1),aa(2),aa(3),xx(1),xx(2),xx(3))
2950 c back to single precision
2951 x1=sngl(xx(1))
2952 x2=sngl(xx(2))
2953 x3=sngl(xx(3))
2954 return
2955 end
2956
2957 c-----------------------------------------------------------------------
2958 subroutine utroa2(phi,a1,a2,a3,x1,x2,x3)
2959 c-----------------------------------------------------------------------
2960 c rotates x by angle phi around axis a.
2961 c normalization of a is irrelevant.
2962 c double precision phi,a1,a2,a3,x1,x2,x3
2963 c-----------------------------------------------------------------------
2964 double precision phi,a1,a2,a3,x1,x2,x3
2965 double precision aaa,aa(3),xxx,xx(3),e1(3),e2(3),e3(3),xp,xt,dphi
2966 dphi=phi
2967 xx(1)=x1
2968 xx(2)=x2
2969 xx(3)=x3
2970 aa(1)=a1
2971 aa(2)=a2
2972 aa(3)=a3
2973 aaa=0d0
2974 xxx=0d0
2975 do i=1,3
2976 aaa=aaa+aa(i)**2
2977 xxx=xxx+xx(i)**2
2978 enddo
2979 if(xxx.eq.0d0)return
2980 if(aaa.eq.0d0)call utstop('utroa1: zero rotation axis&')
2981 aaa=1.0/dsqrt(aaa)
2982 c e3 = a / !a!
2983 do i=1,3
2984 e3(i)=aa(i)*aaa
2985 enddo
2986 c x_parallel
2987 xp=0
2988 do i=1,3
2989 xp=xp+xx(i)*e3(i)
2990 enddo
2991 c x_transverse
2992 if(xxx-xp**2.le.0.)return
2993 xt=dsqrt(xxx-xp**2)
2994 c e1 = vector x_transverse / absolute value x_transverse
2995 do i=1,3
2996 e1(i)=(xx(i)-e3(i)*xp)/xt
2997 enddo
2998 c e2 orthogonal e3,e1
2999 call utvec2(e3,e1,e2)
3000 c rotate x
3001 do i=1,3
3002 xx(i)=xp*e3(i)+xt*cos(dphi)*e1(i)+xt*sin(dphi)*e2(i)
3003 enddo
3004 xxx=0d0
3005 do i=1,3
3006 xxx=xxx+xx(i)**2
3007 enddo
3008 c back to single precision
3009 x1=xx(1)
3010 x2=xx(2)
3011 x3=xx(3)
3012 return
3013 end
3014
3015 cc--------------------------------------------------------------------
3016 c function utroot(funcd,x1,x2,xacc)
3017 cc--------------------------------------------------------------------
3018 cc combination of newton-raphson and bisection method for root finding
3019 cc input:
3020 cc funcd: subr returning fctn value and first derivative
3021 cc x1,x2: x-interval
3022 cc xacc: accuracy
3023 cc output:
3024 cc utroot: root
3025 cc--------------------------------------------------------------------
3026 c include 'epos.inc'
3027 c INTEGER MAXIT
3028 c REAL utroot,x1,x2,xacc
3029 c EXTERNAL funcd
3030 c PARAMETER (MAXIT=100)
3031 c INTEGER j
3032 c REAL df,dx,dxold,f,fh,fl,temp,xh,xl
3033 c call funcd(x1,fl,df)
3034 c call funcd(x2,fh,df)
3035 c if((fl.gt.0..and.fh.gt.0.).or.(fl.lt.0..and.fh.lt.0.))
3036 c *call utstop('utroot: root must be bracketed&')
3037 c if(fl.eq.0.)then
3038 c utroot=x1
3039 c return
3040 c else if(fh.eq.0.)then
3041 c utroot=x2
3042 c return
3043 c else if(fl.lt.0.)then
3044 c xl=x1
3045 c xh=x2
3046 c else
3047 c xh=x1
3048 c xl=x2
3049 c endif
3050 c utroot=.5*(x1+x2)
3051 c dxold=abs(x2-x1)
3052 c dx=dxold
3053 c call funcd(utroot,f,df)
3054 c do 11 j=1,MAXIT
3055 c if(((utroot-xh)*df-f)*((utroot-xl)*df-f).ge.0..or. abs(2.*
3056 c *f).gt.abs(dxold*df) ) then
3057 c dxold=dx
3058 c dx=0.5*(xh-xl)
3059 c utroot=xl+dx
3060 c if(xl.eq.utroot)return
3061 c else
3062 c dxold=dx
3063 c dx=f/df
3064 c temp=utroot
3065 c utroot=utroot-dx
3066 c if(temp.eq.utroot)return
3067 c endif
3068 c if(abs(dx).lt.xacc) return
3069 c call funcd(utroot,f,df)
3070 c if(f.lt.0.) then
3071 c xl=utroot
3072 c else
3073 c xh=utroot
3074 c endif
3075 c11 continue
3076 c call utmsg('utroot')
3077 c write(ifch,*)'***** exceeding maximum iterations'
3078 c write(ifch,*)'dx:',dx
3079 c call utmsgf
3080 c return
3081 c END
3082 c
3083 c-----------------------------------------------------------------------
3084 subroutine utrot2(isig,ax,ay,az,x,y,z)
3085 c-----------------------------------------------------------------------
3086 c performs a rotation, double prec.
3087 c-----------------------------------------------------------------------
3088 include 'epos.inc'
3089 double precision ax,ay,az,x,y,z,rx,ry,rz
3090 *,alp,bet,cosa,sina,cosb,sinb,xs,ys,zs
3091 if(ax**2.eq.0.and.ay**2.eq.0.and.az**2.eq.0.)then
3092 write(ifch,*)'ax**2,ay**2,az**2:',ax**2,ay**2,az**2
3093 write(ifch,*)'ax,ay,az:',ax,ay,az
3094 call utstop('utrot2: zero vector.&')
3095 endif
3096 if(az.ge.0.)then
3097 rx=ax
3098 ry=ay
3099 rz=az
3100 else
3101 rx=-ax
3102 ry=-ay
3103 rz=-az
3104 endif
3105 if(rz**2+ry**2.ne.0.)then
3106 alp=dabs(dacos(rz/dsqrt(rz**2+ry**2)))*sign(1.,sngl(ry))
3107 bet=
3108 *dabs(dacos(dsqrt(rz**2+ry**2)/dsqrt(rz**2+ry**2+rx**2)))*
3109 *sign(1.,sngl(rx))
3110 else
3111 alp=3.1415927d0/2d0
3112 bet=3.1415927d0/2d0
3113 endif
3114 cosa=dcos(alp)
3115 sina=dsin(alp)
3116 cosb=dcos(bet)
3117 sinb=dsin(bet)
3118 if(isig.ge.0)then
3119 xs=x*cosb-y*sina*sinb-z*cosa*sinb
3120 ys= y*cosa -z*sina
3121 zs=x*sinb+y*sina*cosb+z*cosa*cosb
3122 else !if(isig.lt.0)then
3123 xs= x*cosb +z*sinb
3124 ys=-x*sinb*sina+y*cosa+z*cosb*sina
3125 zs=-x*sinb*cosa-y*sina+z*cosb*cosa
3126 endif
3127 x=xs
3128 y=ys
3129 z=zs
3130 return
3131 end
3132
3133
3134 c-----------------------------------------------------------------------
3135 subroutine utrot4(isig,ax,ay,az,x,y,z)
3136 c-----------------------------------------------------------------------
3137 c performs a rotation, double prec.
3138 c arguments partly single
3139 c-----------------------------------------------------------------------
3140 double precision ax,ay,az,xx,yy,zz
3141 xx=dble(x)
3142 yy=dble(y)
3143 zz=dble(z)
3144 call utrot2(isig,ax,ay,az,xx,yy,zz)
3145 x=sngl(xx)
3146 y=sngl(yy)
3147 z=sngl(zz)
3148 return
3149 end
3150
3151 c-----------------------------------------------------------------------
3152 subroutine utrota(isig,ax,ay,az,x,y,z)
3153 c-----------------------------------------------------------------------
3154 c performs a rotation
3155 c-----------------------------------------------------------------------
3156 if(az.ge.0.)then
3157 rx=ax
3158 ry=ay
3159 rz=az
3160 else
3161 rx=-ax
3162 ry=-ay
3163 rz=-az
3164 endif
3165 if(rz.eq.0..and.ry.eq.0.)then
3166 alp=0.
3167 stop
3168 else
3169 alp=abs(utacos(rz/sqrt(rz**2+ry**2)))*sign(1.,ry)
3170 endif
3171 bet=
3172 *abs(utacos(sqrt(rz**2+ry**2)/sqrt(rz**2+ry**2+rx**2)))*sign(1.,rx)
3173 cosa=cos(alp)
3174 sina=sin(alp)
3175 cosb=cos(bet)
3176 sinb=sin(bet)
3177 if(isig.ge.0)then
3178 xs=x*cosb-y*sina*sinb-z*cosa*sinb
3179 ys= y*cosa -z*sina
3180 zs=x*sinb+y*sina*cosb+z*cosa*cosb
3181 else !if(isig.lt.0)then
3182 xs= x*cosb +z*sinb
3183 ys=-x*sinb*sina+y*cosa+z*cosb*sina
3184 zs=-x*sinb*cosa-y*sina+z*cosb*cosa
3185 endif
3186 x=xs
3187 y=ys
3188 z=zs
3189 return
3190 end
3191
3192 c-----------------------------------------------------------------------
3193 subroutine utstop(text)
3194 c-----------------------------------------------------------------------
3195 c returns error message and stops execution.
3196 c text is an optonal text to appear in the error message.
3197 c text is a character string of length 40;
3198 c for shorter text, it has to be terminated by &;
3199 c example: call utstop('error in subr xyz&')
3200 c-----------------------------------------------------------------------
3201 include 'epos.inc'
3202 c parameter(itext=40)
3203 character text*(*) ,txt*6
3204 imax=index(text,'&')
3205 do 1 j=1,2
3206 if(j.eq.1)then
3207 ifi=ifch
3208 else !if(j.eq.2)
3209 ifi=ifmt
3210 endif
3211 if(imax.gt.1)then
3212 write(ifi,101)('*',k=1,72),text(1:imax-1)
3213 *,nrevt+1,nint(seedj),seedc,('*',k=1,72)
3214 else
3215 write(ifi,101)('*',k=1,72),' '
3216 *,nrevt+1,nint(seedj),seedc,('*',k=1,72)
3217 endif
3218 101 format(
3219 *1x,72a1
3220 */1x,'***** stop in ',a
3221 */1x,'***** current event number: ',i12
3222 */1x,'***** initial seed for current run:',i10
3223 */1x,'***** initial seed for current event:',d25.15
3224 */1x,72a1)
3225 1 continue
3226 c c=0.
3227 c b=a/c
3228 stop
3229 entry utmsg(txt)
3230 imsg=imsg+1
3231 write(ifch,'(1x,74a1)')('*',j=1,72)
3232 write(ifch,100)txt,nrevt+1,nint(seedj),seedc
3233 100 format(1x,'***** msg from ',a6,'. es:',i7,2x,i9,2x,d23.17)
3234 return
3235 entry utmsgf
3236 if(ish.eq.1)return
3237 write(ifch,'(1x,74a1)')('*',j=1,72)
3238 end
3239
3240 c-----------------------------------------------------------------
3241 subroutine uttrap(func,a,b,s)
3242 c-----------------------------------------------------------------
3243 c trapezoidal method for integration.
3244 c input: fctn func and limits a,b
3245 c output: value s of the integral
3246 c-----------------------------------------------------------------
3247 include 'epos.inc'
3248
3249 INTEGER JMAX
3250 REAL a,b,func,s
3251 EXTERNAL func
3252 PARAMETER (JMAX=10)
3253 CU USES uttras
3254 INTEGER j
3255 REAL olds
3256 olds=-1.e30
3257 do 11 j=1,JMAX
3258 if(ish.ge.9)write(ifch,*)'sr uttrap: j:',j
3259 call uttras(func,a,b,s,j)
3260 ds=abs(s-olds)
3261 if (ds.lt.epsr*abs(olds)) return
3262 olds=s
3263 11 continue
3264
3265 c-c nepsr=nepsr+1
3266 if(ish.ge.9)then
3267 call utmsg('uttrap')
3268 write(ifch,*)
3269 *'***** requested accuracy could not be achieved'
3270 write(ifch,*)'achieved accuracy: ',ds/abs(olds)
3271 write(ifch,*)'requested accuracy:',epsr
3272 call utmsgf
3273 endif
3274
3275 END
3276
3277 c-----------------------------------------------------------------
3278 subroutine uttraq(func,a,b,s)
3279 c-----------------------------------------------------------------
3280 c trapezoidal method for integration.
3281 c input: function func and limits a,b
3282 c output: value s of the integral
3283 c-----------------------------------------------------------------
3284
3285 REAL a,b,func,s
3286 EXTERNAL func
3287 PARAMETER (eps=1.e-6)
3288 CU USES uttras
3289 INTEGER j
3290 REAL olds
3291 olds=-1.e30
3292 j=1
3293 10 call uttras(func,a,b,s,j)
3294 ds=abs(s-olds)
3295 if (ds.le.eps*abs(olds)) return
3296 olds=s
3297 if(j.ge.15)return
3298 j=j+1
3299 goto10
3300 END
3301
3302 c-----------------------------------------------------------------
3303 subroutine uttras(func,a,b,s,n)
3304 c-----------------------------------------------------------------
3305 c performs one iteration of the trapezoidal method for integration
3306 c-----------------------------------------------------------------
3307 INTEGER n
3308 REAL a,b,s,func
3309 EXTERNAL func
3310 INTEGER it,j
3311 REAL del,sum,tnm,x
3312 if (n.eq.1) then
3313 s=0.5*(b-a)*(func(a)+func(b))
3314 else
3315 it=2**(n-2)
3316 tnm=it
3317 del=(b-a)/tnm
3318 x=a+0.5*del
3319 sum=0.
3320 do 11 j=1,it
3321 sum=sum+func(x)
3322 x=x+del
3323 11 continue
3324 s=0.5*(s+(b-a)*sum/tnm)
3325 endif
3326 return
3327 END
3328
3329 cc-----------------------------------------------------------------------
3330 c subroutine utvec1(a,b,c)
3331 cc-----------------------------------------------------------------------
3332 cc returns vector product c = a x b .
3333 cc-----------------------------------------------------------------------
3334 c real a(3),b(3),c(3)
3335 c c(1)=a(2)*b(3)-a(3)*b(2)
3336 c c(2)=a(3)*b(1)-a(1)*b(3)
3337 c c(3)=a(1)*b(2)-a(2)*b(1)
3338 c return
3339 c end
3340 c
3341 c-----------------------------------------------------------------------
3342 subroutine utvec2(a,b,c)
3343 c-----------------------------------------------------------------------
3344 c returns vector product c = a x b .
3345 c a,b,c double precision.
3346 c-----------------------------------------------------------------------
3347 double precision a(3),b(3),c(3)
3348 c(1)=a(2)*b(3)-a(3)*b(2)
3349 c(2)=a(3)*b(1)-a(1)*b(3)
3350 c(3)=a(1)*b(2)-a(2)*b(1)
3351 return
3352 end
3353
3354 c-------------------------------------------------------------------
3355 subroutine utword(line,i,j,iqu)
3356 c-------------------------------------------------------------------
3357 c finds the first word of the character string line(j+1:1000).
3358 c the word is line(i:j) (with new i and j).
3359 c if j<0 or if no word found --> new line read.
3360 c a text between quotes "..." is considered one word;
3361 c stronger: a text between double quotes ""..."" is consid one word
3362 c stronger: a text between "{ and }" is considered one word
3363 c-------------------------------------------------------------------
3364 c input:
3365 c line: character string (*1000)
3366 c i: integer between 1 and 1000
3367 c iqu: for iqu=1 a "?" is written to output before reading a line,
3368 c otherwise (iqu/=1) nothing is typed
3369 c output:
3370 c i,j: left and right end of word (word=line(i:j))
3371 c-------------------------------------------------------------------
3372 include 'epos.inc'
3373 parameter(mempty=2)
3374 character*1 empty(mempty),mk
3375 character line*1000
3376 character*2 mrk
3377 data empty/' ',','/
3378 parameter(mxdefine=40)
3379 character w1define*100,w2define*100
3380 common/cdefine/ndefine,l1define(mxdefine),l2define(mxdefine)
3381 & ,w1define(mxdefine),w2define(mxdefine)
3382
3383 j0=0
3384 if(j.ge.0)then
3385 i=j
3386 goto 1
3387 endif
3388
3389 5 continue
3390 if(iqu.eq.1.and.iprmpt.gt.0)write(ifmt,'(a)')'?'
3391 if(nopen.eq.0)then
3392 ifopx=ifop
3393 elseif(nopen.gt.0)then
3394 ifopx=20+nopen
3395 else !if(nopen.lt.0)
3396 ifopx=ifcp
3397 endif
3398 read(ifopx,'(a1000)',end=9999)line
3399 if(iecho.eq.1.or.(nopen.ge.0.and.kcpopen.eq.1))then
3400 kmax=2
3401 do k=3,1000
3402 if(line(k:k).ne.' ')kmax=k
3403 enddo
3404 else
3405 kmax=0
3406 endif
3407 if(nopen.ge.0.and.kcpopen.eq.1)
3408 & write(ifcp,'(a)')line(1:kmax)
3409 if(iecho.eq.1)
3410 & write(ifmt,'(a)')line(1:kmax)
3411 i=0
3412
3413 1 i=i+1
3414 if(i.gt.1000)goto 5
3415 if(line(i:i).eq.'!')goto 5
3416 do ne=1,mempty
3417 if(line(i:i).eq.empty(ne))goto 1
3418 enddo
3419
3420 nbla=1
3421 mrk=' '
3422 mk=' '
3423 if(line(i:i).eq.'~')mk='~'
3424 if(line(i:i+1).eq.'"{')mrk='}"'
3425 if(line(i:i+1).eq.'""')mrk='""'
3426 if(mrk.ne.' ')goto 10
3427 if(line(i:i).eq.'"')mk='"'
3428 if(mk.ne.' ')goto 8
3429 j=i-1
3430 6 j=j+1
3431 if(j.gt.1000)goto 7
3432 if(line(j:j).eq.'!')goto 7
3433 do ne=1,mempty
3434 if(line(j:j).eq.empty(ne))goto 7
3435 enddo
3436 goto 6
3437
3438 8 continue
3439 if(i.ge.1000-1)stop'utword: make line shorter!!! '
3440 i=i+1
3441 j=i
3442 if(line(j:j).eq.mk)stop'utword: empty string!!! '
3443 9 j=j+1
3444 if(j.gt.1000)then !reach the end of the line
3445 j=j-nbla+2
3446 goto 7
3447 endif
3448 if(line(j:j).eq.' ')then
3449 nbla=nbla+1
3450 else
3451 nbla=2
3452 endif
3453 if(line(j:j).eq.mk)then
3454 line(i-1:i-1)=' '
3455 line(j:j)=' '
3456 goto 7
3457 endif
3458 goto 9
3459
3460 10 continue
3461 if(i.ge.1000-3)stop'utword: make line shorter!!!! '
3462 i=i+2
3463 j=i
3464 if(line(j:j+1).eq.mrk)stop'utword: empty string!!!! '
3465 11 j=j+1
3466 if(j.gt.1000-1)then !reach the end of the line
3467 j=j-nbla+2
3468 goto 7
3469 endif
3470 if(line(j:j+1).eq.mrk)then
3471 line(i-2:i-1)=' '
3472 line(j:j+1)=' '
3473 goto 7
3474 endif
3475 if(line(j:j).eq.' ')then
3476 nbla=nbla+1
3477 else
3478 nbla=2
3479 endif
3480 goto 11
3481
3482 7 j=j-1
3483 !--------#define---------------
3484 if(ndefine.gt.0)then
3485 do ndf=1,ndefine
3486 l1=l1define(ndf)
3487 l2=l2define(ndf)
3488 do i0=i,j+1-l1
3489 if(line(i0:i0-1+l1).eq.w1define(ndf)(1:l1))then
3490 if(l2.eq.l1)then
3491 line(i0:i0-1+l1)=w2define(ndf)(1:l2)
3492 elseif(l2.lt.l1)then
3493 line(i0:i0+l2-1)=w2define(ndf)(1:l2)
3494 do k=i0+l2,i0-1+l1
3495 line(k:k)=' '
3496 enddo
3497 elseif(l2.gt.l1)then
3498 do k=i0+l1,i0+l2-1
3499 if(line(k:k).ne.' ')
3500 & stop'utword: no space for `define` replacement. '
3501 enddo
3502 line(i0:i0+l2-1)=w2define(ndf)(1:l2)
3503 j=i0+l2-1
3504 endif
3505 endif
3506 enddo
3507 enddo
3508 do k=i,j
3509 if(line(k:k).ne.' ')j0=j
3510 enddo
3511 j=j0
3512 endif
3513 !--------
3514 return
3515
3516 9999 close(ifopx)
3517 nopen=nopen-1
3518 if(nopen.eq.0.and.iprmpt.eq.-1)iprmpt=1
3519 goto 5
3520 end
3521
3522 c--------------------------------------------------------------------
3523 subroutine utworn(line,j,ne)
3524 c--------------------------------------------------------------------
3525 c returns number ne of nonempty characters of line(j+1:1000)
3526 c--------------------------------------------------------------------
3527 character line*1000
3528 ne=0
3529 do l=j+1,1000
3530 if(line(l:l).ne.' ')ne=ne+1
3531 enddo
3532 return
3533 end
3534
3535 c-----------------------------------------------------------------------
3536 subroutine getairmol(iz,ia)
3537 c-----------------------------------------------------------------------
3538 include 'epos.inc'
3539 i=0
3540 r=rangen()
3541 do while(r.gt.0.) ! choose air-molecule
3542 i=i+1
3543 r=r-airwnxs(i)
3544 enddo
3545 iz = nint(airznxs(i))
3546 ia = nint(airanxs(i))
3547 end
3548
3549 c----------------------------------------------------------------------
3550
3551 subroutine factoriel
3552
3553 c----------------------------------------------------------------------
3554 c tabulation of fctrl(n)=n!, facto(n)=1/n! and utgamtab(x) for x=0 to 50
3555 c----------------------------------------------------------------------
3556 include 'epos.incems'
3557 double precision utgamtab,utgam,x
3558 common/gamtab/utgamtab(10000)
3559
3560 nfctrl=100
3561 fctrl(0)=1.D0
3562 facto(0)=1.D0
3563 do i=1,min(npommx,nfctrl)
3564 fctrl(i)=fctrl(i-1)*dble(i)
3565 facto(i)=1.d0/fctrl(i)
3566 enddo
3567
3568 do k=1,10000
3569 x=dble(k)/100.d0
3570 utgamtab(k)=utgam(x)
3571 enddo
3572
3573 return
3574 end
3575
3576 c-----------------------------------------------------------------------
3577 subroutine fremnu(amin,ca,cb,ca0,cb0,ic1,ic2,ic3,ic4)
3578 c-----------------------------------------------------------------------
3579 common/hadr2/iomodl,idproj,idtarg,wexcit
3580 real pnll,ptq
3581 common/hadr1/pnll,ptq,exmass,cutmss,wproj,wtarg
3582 parameter (nflav=6)
3583 integer ic(2),jc(nflav,2)
3584 ic(1)=ca
3585 ic(2)=cb
3586 call iddeco(ic,jc)
3587 keu=jc(1,1)-jc(1,2)
3588 ked=jc(2,1)-jc(2,2)
3589 kes=jc(3,1)-jc(3,2)
3590 kec=jc(4,1)-jc(4,2)
3591 keb=jc(5,1)-jc(5,2)
3592 ket=jc(6,1)-jc(6,2)
3593 amin=utamnu(keu,ked,kes,kec,keb,ket,5) !???4=2mults, 5=1mult
3594 if(ca-ca0.eq.0..and.cb-cb0.eq.0..and.rangen().gt.wexcit)then
3595 ic3=0
3596 ic4=0
3597 ic1=ca
3598 ic2=cb
3599 else
3600 amin=amin+exmass
3601 n=0
3602 do i=1,4
3603 do j=1,2
3604 n=n+jc(i,j)
3605 enddo
3606 enddo
3607 k=1+rangen()*n
3608 do i=1,4
3609 do j=1,2
3610 k=k-jc(i,j)
3611 if(k.le.0)goto 1
3612 enddo
3613 enddo
3614 1 if(j.eq.1)then
3615 ic3=10**(6-i)
3616 ic4=0
3617 else
3618 ic3=0
3619 ic4=10**(6-i)
3620 endif
3621 ic1=int(ca)-ic3
3622 ic2=int(cb)-ic4
3623 endif
3624 return
3625 end
3626
3627
3628 c-----------------------------------------------------------------------
3629 function fremnux(jc)
3630 c-----------------------------------------------------------------------
3631 real pnll,ptq
3632 common/hadr1/pnll,ptq,exmass,cutmss,wproj,wtarg
3633 parameter (nflav=6)
3634 integer jc(nflav,2)!,ic(2)
3635 c ic(1)=210000
3636 c ic(2)=0
3637 c call iddeco(ic,jc)
3638 keu=jc(1,1)-jc(1,2)
3639 ked=jc(2,1)-jc(2,2)
3640 kes=jc(3,1)-jc(3,2)
3641 kec=jc(4,1)-jc(4,2)
3642 keb=jc(5,1)-jc(5,2)
3643 ket=jc(6,1)-jc(6,2)
3644 fremnux=utamnu(keu,ked,kes,kec,keb,ket,4) !+exmass !???4=2mults, 5=1mult
3645 return
3646 end
3647
3648 c-----------------------------------------------------------------------
3649 function fremnux2(jc)
3650 c-----------------------------------------------------------------------
3651 real pnll,ptq
3652 common/hadr1/pnll,ptq,exmass,cutmss,wproj,wtarg
3653 parameter (nflav=6)
3654 integer jc(nflav,2)!,ic(2)
3655 c ic(1)=210000
3656 c ic(2)=0
3657 c call iddeco(ic,jc)
3658 keu=jc(1,1)-jc(1,2)
3659 ked=jc(2,1)-jc(2,2)
3660 kes=jc(3,1)-jc(3,2)
3661 kec=jc(4,1)-jc(4,2)
3662 keb=jc(5,1)-jc(5,2)
3663 ket=jc(6,1)-jc(6,2)
3664 fremnux2=utamnu(keu,ked,kes,kec,keb,ket,5) !+exmass !???4=2mults, 5=1mult
3665 return
3666 end
3667
3668 c-----------------------------------------------------------------------
3669 function fremnux3(jci)
3670 c-----------------------------------------------------------------------
3671 c minimum mass from ic counting all quarks
3672 c-----------------------------------------------------------------------
3673 include 'epos.inc'
3674 integer jc(nflav,2),jci(nflav,2)!,ic(2)
3675 c ic(1)=210000
3676 c ic(2)=0
3677 c print *,'start',ic
3678 fremnux3=0.
3679 do j=1,2
3680 do i=1,nflav
3681 jc(i,j)=jci(i,j)
3682 enddo
3683 enddo
3684 c call iddeco(ic,jc)
3685 call idquacjc(jc,nqua,naqu)
3686 do ii=1,2
3687 if(ii.eq.1)then
3688 nqu=nqua
3689 else
3690 nqu=naqu
3691 endif
3692 if(nqu.ge.3)then
3693 do while(jc(3,ii).ne.0.and.nqu.ge.3) !count baryons with s quark
3694 jc(3,ii)=jc(3,ii)-1
3695 if(jc(3,ii).gt.0)then
3696 jc(3,ii)=jc(3,ii)-1
3697 if(jc(3,ii).gt.0)then
3698 jc(3,ii)=jc(3,ii)-1
3699 fremnux3=fremnux3+asuhax(4)
3700 elseif(jc(2,ii).gt.0)then
3701 jc(2,ii)=jc(2,ii)-1
3702 fremnux3=fremnux3+asuhax(4)
3703 elseif(jc(1,ii).gt.0)then
3704 jc(1,ii)=jc(1,ii)-1
3705 fremnux3=fremnux3+asuhax(4)
3706 endif
3707 elseif(jc(2,ii).gt.0)then
3708 jc(2,ii)=jc(2,ii)-1
3709 if(jc(1,ii).gt.0)then
3710 jc(1,ii)=jc(1,ii)-1
3711 fremnux3=fremnux3+asuhax(3)
3712 elseif(jc(2,ii).gt.0)then
3713 jc(2,ii)=jc(2,ii)-1
3714 fremnux3=fremnux3+asuhax(3)
3715 endif
3716 elseif(jc(1,ii).gt.0)then
3717 jc(1,ii)=jc(1,ii)-2
3718 fremnux3=fremnux3+asuhay(3)
3719 endif
3720 nqu=nqu-3
3721 enddo
3722 do while(jc(2,ii).ne.0.and.nqu.ge.3) !count baryons with d quark
3723 jc(2,ii)=jc(2,ii)-1
3724 if(jc(1,ii).gt.0)then
3725 jc(1,ii)=jc(1,ii)-1
3726 if(jc(2,ii).gt.0)then
3727 jc(2,ii)=jc(2,ii)-1
3728 fremnux3=fremnux3+asuhay(2)
3729 elseif(jc(1,ii).gt.0)then
3730 jc(1,ii)=jc(1,ii)-1
3731 fremnux3=fremnux3+asuhay(2)
3732 endif
3733 elseif(jc(2,ii).gt.0)then
3734 jc(2,ii)=jc(2,ii)-2
3735 fremnux3=fremnux3+asuhay(3)
3736 endif
3737 nqu=nqu-3
3738 enddo
3739 do while(jc(1,ii).ne.0.and.nqu.ge.3) !count baryons with s quark
3740 jc(1,ii)=jc(1,ii)-3
3741 fremnux3=fremnux3+asuhay(3)
3742 nqu=nqu-3
3743 enddo
3744 if(ii.eq.1)then
3745 nqua=nqu
3746 else
3747 naqu=nqu
3748 endif
3749 endif
3750 c print *,ii,nqua,naqu,jc,fremnux3
3751 enddo
3752 if(nqua+naqu.ne.0)then
3753 do while(jc(3,1).ne.0) !count mesons with s quark
3754 jc(3,1)=jc(3,1)-1
3755 if(jc(3,2).gt.0)then
3756 jc(3,2)=jc(3,2)-1
3757 fremnux3=fremnux3+asuhax(6)
3758 elseif(jc(2,2).gt.0)then
3759 jc(2,2)=jc(2,2)-1
3760 fremnux3=fremnux3+asuhay(6)
3761 elseif(jc(1,2).gt.0)then
3762 jc(1,2)=jc(1,2)-1
3763 fremnux3=fremnux3+asuhay(6)
3764 endif
3765 enddo
3766 do while(jc(2,1).ne.0) !count mesons with d quark
3767 jc(2,1)=jc(2,1)-1
3768 if(jc(2,2).gt.0)then
3769 jc(2,2)=jc(2,2)-1
3770 fremnux3=fremnux3+asuhay(5)
3771 elseif(jc(1,2).gt.0)then
3772 jc(1,2)=jc(1,2)-1
3773 fremnux3=fremnux3+asuhay(5)
3774 endif
3775 enddo
3776 do while(jc(1,1).ne.0) !count mesons with s quark
3777 jc(1,1)=jc(1,1)-1
3778 if(jc(1,2).gt.0)then
3779 jc(1,2)=jc(1,2)-1
3780 fremnux3=fremnux3+asuhay(5)
3781 endif
3782 enddo
3783 endif
3784 c fremnux3=fremnux3+0.5
3785 c print *,'stop',nqua,naqu,fremnux3
3786
3787 return
3788 end
3789
3790 c-----------------------------------------------------------------------
3791 subroutine fremnx(ammax,amin,sm,ic3,ic4,iret)
3792 c-----------------------------------------------------------------------
3793 common/psar9/ alpr
3794 include 'epos.inc'
3795 iret=0
3796 if(ic3.eq.0.and.ic4.eq.0)then
3797 if(ammax.lt.amin**2)then
3798 iret=1
3799 return
3800 endif
3801 sm=amin**2
3802 else
3803 c ammax1=min(ammax,(engy/4.)**2)
3804 ammax1=ammax
3805 if(ammax1.lt.amin**2)then
3806 iret=1
3807 return
3808 endif
3809 if(alpr.eq.-1.)then
3810 sm=amin**2*(ammax1/amin**2)**rangen()
3811 else
3812 sm=amin**2*(1.+((ammax1/amin**2)**(1.+alpr)-1.)
3813 * *rangen())**(1./(1.+alpr))
3814 endif
3815 endif
3816 return
3817 end
3818
3819 SUBROUTINE gaulag(x,w,n,alf)
3820 INTEGER n,MAXIT
3821 REAL alf,w(n),x(n)
3822 DOUBLE PRECISION EPS
3823 PARAMETER (EPS=3.D-14,MAXIT=10)
3824 CU USES gammln
3825 INTEGER i,its,j
3826 REAL ai,gammln
3827 DOUBLE PRECISION p1,p2,p3,pp,z,z1
3828 z=0.
3829 do 13 i=1,n
3830 if(i.eq.1)then
3831 z=(1.+alf)*(3.+.92*alf)/(1.+2.4*n+1.8*alf)
3832 else if(i.eq.2)then
3833 z=z+(15.+6.25*alf)/(1.+.9*alf+2.5*n)
3834 else
3835 ai=i-2
3836 z=z+((1.+2.55*ai)/(1.9*ai)+1.26*ai*alf/(1.+3.5*ai))*
3837 *(z-x(i-2))/(1.+.3*alf)
3838 endif
3839 do 12 its=1,MAXIT
3840 p1=1.d0
3841 p2=0.d0
3842 do 11 j=1,n
3843 p3=p2
3844 p2=p1
3845 p1=((2*j-1+alf-z)*p2-(j-1+alf)*p3)/j
3846 11 continue
3847 pp=(n*p1-(n+alf)*p2)/z
3848 z1=z
3849 z=z1-p1/pp
3850 if(abs(z-z1).le.EPS)goto 1
3851 12 continue
3852 call utstop("too many iterations in gaulag")
3853 1 x(i)=z
3854 w(i)=-exp(gammln(alf+n)-gammln(float(n)))/(pp*n*p2)
3855 13 continue
3856 return
3857 END
3858 C (C) Copr. 1986-92 Numerical Recipes Software 4+1$!].
3859
3860 FUNCTION gammln(xx)
3861 REAL gammln,xx
3862 INTEGER j
3863 DOUBLE PRECISION ser,stp,tmp,x,y,cof(6)
3864 SAVE cof,stp
3865 DATA cof,stp/76.18009172947146d0,-86.50532032941677d0,
3866 *24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,
3867 *-.5395239384953d-5,2.5066282746310005d0/
3868 x=xx
3869 y=x
3870 tmp=x+5.5d0
3871 tmp=(x+0.5d0)*log(tmp)-tmp
3872 ser=1.000000000190015d0
3873 do 11 j=1,6
3874 y=y+1.d0
3875 ser=ser+cof(j)/y
3876 11 continue
3877 gammln=tmp+log(stp*ser/x)
3878 return
3879 END
3880 C (C) Copr. 1986-92 Numerical Recipes Software 4+1$!].
3881
3882 function polar(x,y)
3883
3884 pi=3.1415927
3885 if(abs(x).gt.1.e-6)then
3886 phi=atan(y/x)
3887 if(x.lt.0.)phi=pi+phi
3888 if(phi.lt.0)phi=2*pi+phi
3889 else
3890 phi=0.5*pi
3891 if(y.lt.0)phi=phi+pi
3892 endif
3893 polar=phi
3894
3895 end
3896
3897 subroutine getJKNcentr
3898 end
