File indexing completed on 2024-04-06 12:14:00
0001
0002 subroutine conaa(iret)
0003
0004
0005
0006
0007 include 'epos.inc'
0008 include 'epos.incems'
0009 include 'epos.incsem'
0010 include 'epos.incpar'
0011
0012 common/geom/rmproj,rmtarg,bmax,bkmx
0013 common/nucl3/phi,bimp
0014 common /cncl/xproj(mamx),yproj(mamx),zproj(mamx)
0015 * ,xtarg(mamx),ytarg(mamx),ztarg(mamx)
0016 common/cfacmss/facmss
0017 double precision yyrmax
0018 common/scrangle/ phik3(kollmx),thetak3(kollmx)
0019
0020 call utpri('conaa ',ish,ishini,4)
0021
0022 iret=0
0023
0024
0025
0026
0027 vel=tanh(ypjtl-yhaha)+tanh(yhaha)
0028
0029
0030
0031
0032 if(iokoll.eq.1)then
0033
0034 koll=matarg
0035 do k=1,koll
0036 do n=1,4
0037 coord(n,k)=0.
0038 enddo
0039 enddo
0040 bimp=0
0041 phi=0
0042 xproj(1)=0
0043 yproj(1)=0
0044 zproj(1)=0
0045 lproj(1)=koll
0046 lproj3(1)=0
0047 do k=1,koll
0048 bij=bkmx*sqrt(rangen())
0049 bk(k)=bij
0050 iproj(k)=1
0051 itarg(k)=k
0052 phi=2.*pi*rangen()
0053 xtarg(k)=bij*cos(phi)
0054 ytarg(k)=bij*sin(phi)
0055 ztarg(k)=0
0056 ltarg(k)=1
0057 kproj(1,k)=k
0058 ktarg(k,1)=k
0059 ltarg3(k)=0
0060 ktarg3(k,1)=0
0061 kproj3(k,1)=0
0062 if(iscreen.ne.0.and.bij.le.bkmxndif)then
0063 if(zbrmax.le.0..or. bij.lt.zbcutx+zbrmax*rangen())then
0064 lproj3(1)=lproj3(1)+1
0065 ltarg3(k)=1
0066 kproj3(1,lproj3(1))=k
0067 ktarg3(k,1)=k
0068 endif
0069 endif
0070 enddo
0071
0072 elseif(maproj.eq.1.and.matarg.eq.1)then
0073
0074 b1=bminim
0075 b2=amin1(bkmx,bmaxim)
0076 if(b1.gt.b2)call utstop('conaa: bmin > bmax&')
0077 bimp=sqrt(b1*b1+(b2*b2-b1*b1)*rangen())
0078 koll=1
0079 do n=1,4
0080 coord(n,1)=0.
0081 enddo
0082 bk(1)=bimp
0083 iproj(1)=1
0084 itarg(1)=1
0085 phi=2.*pi*rangen()
0086 xproj(1)=bimp*cos(phi)
0087 yproj(1)=bimp*sin(phi)
0088 zproj(1)=0
0089 xtarg(1)=0
0090 ytarg(1)=0
0091 ztarg(1)=0
0092 lproj(1)=1
0093 ltarg(1)=1
0094 lproj3(1)=1
0095 ltarg3(1)=1
0096 kproj3(1,1)=1
0097 ktarg3(1,1)=1
0098 kproj(1,1)=1
0099 ktarg(1,1)=1
0100
0101 else
0102
0103 call conxyz('p',mamx,xproj,yproj,zproj,ypjtl-yhaha)
0104 call conxyz('t',mamx,xtarg,ytarg,ztarg,yhaha)
0105
0106 bx=0
0107 by=0
0108 if(maproj.gt.0)then
0109 if(bimevt.lt.0)then
0110 b1=bminim
0111 b2=amin1(rmproj+rmtarg,bmaxim)
0112 if(b1.gt.b2)call utstop('conaa: bmin > bmax&')
0113 bimp=sqrt(b1**2+(b2**2-b1**2)*rangen())
0114 if(nbarray.gt.0)bimp=barray(mod(nrevt,nbarray)+1)
0115 if(jpsi.gt.0)then
0116 bimp=b1+(b2-b1)*(float(mod(nrevt,12))+rangen())/12.
0117 bimevt=bimp
0118 endif
0119 phi=phimin+rangen()*(phimax-phimin)
0120 else
0121 phi=phievt
0122 bimp=bimevt
0123 endif
0124 bx=cos(phi)*bimp
0125 by=sin(phi)*bimp
0126 endif
0127 if(jpsi.lt.0)then !modify b
0128 bx=xtarg(1)
0129 by=ytarg(1)
0130 endif
0131 if(maproj.eq.0)goto1000
0132 koll=0
0133 do i=1,maproj
0134 lproj(i)=0
0135 lproj3(i)=0
0136 enddo
0137 do j=1,matarg
0138 ltarg(j)=0
0139 ltarg3(j)=0
0140 enddo
0141 do 12 i=1,maproj
0142 do 11 j=1,matarg
0143 if(jpsi.lt.0.and.ztarg(j).le.ztarg(1))goto11
0144 bij=sqrt((xproj(i)+bx-xtarg(j))**2+(yproj(i)+by-ytarg(j))**2)
0145 if(ish.ge.7)write(ifch,*)'i_p:',i,' i_t:',j,' b_ij:',bij
0146 if(bij.gt.bkmx)goto 11
0147
0148 koll=koll+1
0149 if(koll.gt.kollmx)call utstop('conaa: kollmx too small&')
0150 bk(koll)=bij
0151 bkx(koll)=xproj(i)+bx-xtarg(j)
0152 bky(koll)=yproj(i)+by-ytarg(j)
0153 iproj(koll)=i
0154 itarg(koll)=j
0155 lproj(i)=lproj(i)+1
0156 ltarg(j)=ltarg(j)+1
0157 kproj(i,lproj(i))=koll
0158 ktarg(j,ltarg(j))=koll
0159 phik3(koll)=0.
0160 thetak3(koll)=0.
0161 if(iscreen.ne.0.and.bij.le.bkmxndif)then
0162 if(zbrmax.gt.0..and.bij.gt.zbcutx+zbrmax*rangen())goto 11
0163 lproj3(i)=lproj3(i)+1
0164 ltarg3(j)=ltarg3(j)+1
0165 kproj3(i,lproj3(i))=koll
0166 ktarg3(j,ltarg3(j))=koll
0167
0168 if(abs(bky(koll)).gt.1.e-6)then
0169 if(abs(bkx(koll)).gt.1.e-6)then
0170 phik3(koll)=atan(bky(koll)/bkx(koll))
0171 else
0172 phik3(koll)=sign(0.5*pi,bky(koll))
0173 endif
0174 elseif(bkx(koll).lt.0.)then
0175 phik3(koll)=pi
0176 endif
0177 if(bk(koll).gt.0.)then
0178 thetak3(koll)=atan(bglaubx/bk(koll))
0179 else
0180 thetak3(koll)=0.5*pi
0181 endif
0182 endif
0183
0184 11 continue
0185 12 continue
0186
0187 do k=1,koll
0188 do n=1,4
0189 coord(n,k)=0.
0190 enddo
0191 enddo
0192
0193
0194 endif
0195
0196 if(ish.ge.3)write(ifch,*)'koll=',koll
0197 if(koll.eq.0)goto 1001
0198
0199
0200
0201
0202 do kl=1,koll
0203 i=iproj(kl)
0204 j=itarg(kl)
0205 dist=ztarg(j)-zproj(i)
0206 coord(1,kl)=(xproj(i)+xtarg(j))*0.5
0207 coord(2,kl)=(yproj(i)+ytarg(j))*0.5
0208 coord(3,kl)=(zproj(i)+ztarg(j))*0.5
0209 coord(4,kl)=dist/vel
0210 enddo
0211
0212 if(iscreen.ne.0)call CalcScrPair(bimp)
0213
0214 !~~~~~redefine energy in case of imposed radial flow~~~~~~~~~~~~~~~~
0215 yrmaxi=max(0.,yradmx+yradmi*log(1.+engy/sqrt(float(koll))))
0216 if(yrmaxi.gt.1e-5)then
0217 yyrmax=dble(yrmaxi)
0218 fradflii=sngl(1d0/
0219 & ((sinh(yyrmax)*yyrmax-cosh(yyrmax)+1d0)/(yyrmax**2/2d0)))
0220 else
0221 fradflii=1.
0222 endif
0223 if(ish.ge.3)write(ifch,*)'yrmaxi=',yrmaxi
0224
0225
0226
0227
0228 1000 continue
0229 if(ish.ge.5)then
0230 write(ifch,*)'ia,x/y/zproj:'
0231 do mm=1,maproj
0232 write(ifch,*)mm,xproj(mm),yproj(mm),zproj(mm)
0233 enddo
0234 write(ifch,*)'ia,x/y/ztarg:'
0235 do mm=1,matarg
0236 write(ifch,*)mm,xtarg(mm),ytarg(mm),ztarg(mm)
0237 enddo
0238 write(ifch,*)'iret',iret
0239 endif
0240 call utprix('conaa ',ish,ishini,4)
0241 return
0242
0243 1001 continue !iret=1 causes redo of whole collision
0244 iret=1
0245 if(ish.ge.3)then
0246 write(ifch,*)
0247 write(ifch,*)'***** subroutine conaa:'
0248 write(ifch,*)'***** no nucleon pair found --> no interaction'
0249 write(ifch,*)
0250 endif
0251 goto 1000
0252
0253 end
0254
0255
0256 function frdmzcut(b)
0257
0258
0259
0260
0261 include 'epos.inc'
0262 include 'epos.incpar'
0263
0264 frdmzcut=zbcutx*exp(-(b*rangen()))
0265
0266
0267 return
0268 end
0269
0270
0271 subroutine CalcScrPair(b)
0272
0273
0274
0275
0276 include 'epos.inc'
0277 include 'epos.incems'
0278 include 'epos.incpar'
0279
0280 common/scrangle/ phik3(kollmx),thetak3(kollmx)
0281 logical lascr(kollmx),cont
0282
0283
0284
0285
0286 do i=1,maproj
0287 if(lproj3(i).gt.1)then
0288 do l=2,lproj3(i)
0289 m=l
0290 100 continue
0291 kp=kproj3(i,m-1)
0292 kl=kproj3(i,m)
0293 if(bk(kl).lt.bk(kp))then
0294 kproj3(i,m-1)=kl
0295 kproj3(i,m)=kp
0296 m=m-1
0297 if(m.gt.1)goto 100
0298 endif
0299 enddo
0300 endif
0301 enddo
0302 do j=1,matarg
0303 if(ltarg3(j).gt.1)then
0304 do l=2,ltarg3(j)
0305 m=l
0306 200 continue
0307 kp=ktarg3(j,m-1)
0308 kl=ktarg3(j,m)
0309 if(bk(kl).lt.bk(kp))then
0310 ktarg3(j,m-1)=kl
0311 ktarg3(j,m)=kp
0312 m=m-1
0313 if(m.gt.1)goto 200
0314 endif
0315 enddo
0316 endif
0317 enddo
0318
0319 if(koll.gt.1)then
0320
0321
0322
0323
0324
0325 do i=1,maproj
0326 if(lproj3(i).gt.1)then
0327 kl=kproj3(i,1)
0328 lascr(kl)=.true.
0329 do 300 l=lproj3(i),2,-1
0330 kl=kproj3(i,l)
0331 lascr(kl)=.true.
0332 if(bk(kl).ge.frdmzcut(b))then
0333 do m=1,l-1
0334 km=kproj3(i,m)
0335 if(kl.ne.km.and.bk(km).ge.frdmzcut(b)
0336 & .and.phik3(kl).ge.phik3(km)-thetak3(km)
0337 & .and.phik3(kl).le.phik3(km)+thetak3(km))then
0338 lascr(kl)=.false.
0339 goto 300
0340 endif
0341 enddo
0342 endif
0343 300 continue
0344 endif
0345 enddo
0346
0347 do i=1,maproj
0348 if(lproj3(i).gt.1)then
0349 do l=1,lproj3(i)
0350 kl=kproj3(i,l)
0351 enddo
0352 endif
0353 enddo
0354 do i=1,maproj
0355 if(lproj3(i).gt.1)then
0356 n=2
0357 cont=lproj3(i).gt.1
0358 do while(cont)
0359 l=lproj3(i)
0360 kl=kproj3(i,l)
0361
0362 do while(.not.lascr(kl).and.l.gt.1)
0363 kproj3(i,l)=0
0364 lproj3(i)=lproj3(i)-1
0365 l=lproj3(i)
0366 kl=kproj3(i,l)
0367 enddo
0368 cont=lproj3(i).gt.n
0369 if(cont)then
0370
0371 kn=kproj3(i,n)
0372 if(.not.lascr(kn))then
0373 m=lproj3(i) !last pair always active
0374 km=kproj3(i,m)
0375 kproj3(i,n)=km
0376 kproj3(i,m)=0
0377 lproj3(i)=lproj3(i)-1
0378 endif
0379 endif
0380 n=min(lproj3(i)-1,n)+1
0381 cont=lproj3(i).ne.n
0382 enddo
0383 endif
0384 enddo
0385
0386
0387
0388 do j=1,matarg
0389 if(ltarg3(j).gt.1)then
0390 kl=ktarg3(j,1)
0391 lascr(kl)=.true.
0392 do 400 l=ltarg3(j),2,-1
0393 kl=ktarg3(j,l)
0394 lascr(kl)=.true.
0395 if(bk(kl).ge.frdmzcut(b))then
0396 do m=1,l-1
0397 km=ktarg3(j,m)
0398 if(km.ne.kl.and.bk(km).ge.frdmzcut(b)
0399 & .and.phik3(kl).ge.phik3(km)-thetak3(km)
0400 & .and.phik3(kl).le.phik3(km)+thetak3(km))then
0401 lascr(kl)=.false.
0402 goto 400
0403 endif
0404 enddo
0405 endif
0406 400 continue
0407 endif
0408 enddo
0409
0410 do j=1,matarg
0411 if(ltarg3(j).gt.1)then
0412 n=2
0413 cont=ltarg3(j).gt.1
0414 do while(cont)
0415 l=ltarg3(j)
0416 kl=ktarg3(j,l)
0417
0418 do while(.not.lascr(kl).and.l.gt.1)
0419 ktarg3(j,l)=0
0420 ltarg3(j)=ltarg3(j)-1
0421 l=ltarg3(j)
0422 kl=ktarg3(j,l)
0423 enddo
0424 cont=ltarg3(j).gt.n
0425 if(cont)then
0426
0427 kn=ktarg3(j,n)
0428 if(.not.lascr(kn))then
0429 m=ltarg3(j) !last pair always active
0430 km=ktarg3(j,m)
0431 ktarg3(j,n)=km
0432 ktarg3(j,m)=0
0433 ltarg3(j)=ltarg3(j)-1
0434 endif
0435 endif
0436 n=min(ltarg3(j)-1,n)+1
0437 cont=ltarg3(j).ne.n
0438 enddo
0439 endif
0440 enddo
0441 endif
0442
0443
0444 end
0445
0446
0447 subroutine xGeometry(iii)
0448
0449 include 'epos.inc'
0450 include 'epos.incems'
0451 include 'epos.incsem'
0452 common/xgeom/nnn,naa(kollmx),nbb(kollmx)
0453 character*5 fmt1,fmt2
0454
0455 if(iii.eq.0)then
0456
0457 do k=1,kollmx
0458 naa(k)=0
0459 enddo
0460 nnn=0
0461
0462
0463 elseif(iii.eq.1)then
0464
0465 ngl=0
0466 do k=1,koll
0467 r=bk(k)
0468 if(r.le.sqrt(sigine/10./pi))ngl=ngl+1
0469 enddo
0470 if(ngl.ne.0)then
0471 nnn=nnn+1
0472 naa(ngl)=naa(ngl)+1
0473 endif
0474
0475 elseif(iii.eq.2)then
0476
0477 if(xpar1.eq.0..and.xpar2.eq.0.)then
0478 print*,'---------- xGeometry -----------------------'
0479 return
0480 endif
0481 x1=1-0.01*xpar2
0482 x2=1-0.01*xpar1
0483 kmx=0
0484 nbb(1)=naa(1)
0485 do k=2,kollmx
0486 if(naa(k).ne.0.)kmx=k
0487 nbb(k)=nbb(k-1)+naa(k)
0488 enddo
0489 k1=0
0490 k2=0
0491 do k=1,kmx
0492 x=nbb(k)/float(nnn)
0493 if(x.lt.x1)k1=k
0494 if(x.lt.x2)k2=k
0495 enddo
0496 k1=k1+1
0497 k2=k2+1
0498 x=0
0499 av=0
0500 su=0
0501 p1=0.
0502 p2=0.
0503 do k=1,kmx
0504 xb=x
0505 x=nbb(k)/float(nnn)
0506 y=naa(k)/float(nnn)
0507 dx=x-xb
0508 p=0.
0509 if(k.eq.k1)then
0510 p=(x-x1)/dx
0511 p1=p
0512 elseif(k.eq.k2)then
0513 p=(x2-xb)/dx
0514 p2=p
0515 elseif(k.gt.k1.and.k.lt.k2)then
0516 p=1
0517 endif
0518 av=av+y*p*k
0519 su=su+y*p
0520 enddo
0521 av=av/su
0522 n1=nint(100*p1)
0523 n2=nint(100*p2)
0524 if(n1.eq.0)then
0525 k1=k1+1
0526 n1=100
0527 endif
0528 if(k1.le.9)fmt1='i1,4x'
0529 if(k1.gt.9.and.k1.le.99)fmt1='i2,3x'
0530 if(k1.gt.99.and.k1.le.999)fmt1='i3,2x'
0531 if(k1.gt.999.and.k1.le.9999)fmt1='i4,1x'
0532 if(k2.le.9)fmt2='i1,4x'
0533 if(k2.gt.9.and.k2.le.99)fmt2='i2,3x'
0534 if(k2.gt.99.and.k2.le.999)fmt2='i3,2x'
0535 if(k2.gt.999.and.k2.le.9999)fmt2='i4,1x'
0536 write(6,'(i4,a,i5,a,'//fmt1//',i6,a,i5,a,'//fmt2//',5x,a,f8.2)')
0537 & nint(xpar2),':MIN',n1,'%',k1
0538 & ,nint(xpar1),':MAX',n2,'%',k2 ,'av:',av
0539 endif
0540
0541 end
0542
0543
0544 function conbmx()
0545
0546 double precision om1intbc,p,eps
0547 include 'epos.inc'
0548 include 'epos.incsem'
0549
0550 conbmx=0.
0551 b1=0.
0552 b2=7.
0553 eps=5.0d-3
0554 p=1.d0-dexp(-om1intbc(b2))
0555 if(p.gt.2.d0*eps)return
0556
0557 ntry=0
0558
0559 10 ntry=ntry+1
0560 b=b1+.5*(b2-b1)
0561
0562 p=(1.d0-dexp(-om1intbc(b)))
0563
0564 if(p.gt.eps)then
0565 if(p.gt.2.d0*eps)then
0566 b1=b
0567 else
0568 conbmx=b
0569 return
0570 endif
0571 else
0572 if(p.lt.eps/5.d0)then
0573 b2=b
0574 else
0575 conbmx=b
0576 return
0577 endif
0578 endif
0579
0580 if(ntry.le.1000)goto 10
0581 write(ifmt,*)'Too much try in conbmx ... bmax=',b
0582 conbmx=b
0583 return
0584
0585 end
0586
0587
0588 function conbmxndif()
0589
0590 double precision om1intbc,p,eps
0591 include 'epos.inc'
0592 include 'epos.incsem'
0593
0594
0595 iomegasave=iomega
0596 iomega=2
0597 conbmxndif=0.
0598 b1=0.
0599 b2=7.
0600 conbmxndif=b2
0601 eps=1d-10
0602 p=1.d0-dexp(-om1intbc(b2))
0603 if(p.gt.2.d0*eps)goto 100
0604
0605 ntry=0
0606
0607 10 ntry=ntry+1
0608 b=b1+.5*(b2-b1)
0609
0610 p=(1.d0-dexp(-om1intbc(b)))
0611
0612 if(p.gt.eps)then
0613 if(p.gt.2.d0*eps)then
0614 b1=b
0615 else
0616 conbmxndif=b
0617 goto 100
0618 endif
0619 else
0620 if(p.lt.eps/5.d0)then
0621 b2=b
0622 else
0623 conbmxndif=b
0624 goto 100
0625 endif
0626 endif
0627
0628 if(ntry.le.1000)goto 10
0629 write(ifmt,*)'Too much try in conbmxndif ... bkmxndif=',b
0630 conbmxndif=b
0631 100 iomega=iomegasave
0632 return
0633
0634 end
0635
0636
0637 function conbmxdif()
0638
0639
0640
0641 double precision om1intbc,pmax,drootom,pdiff
0642 include 'epos.inc'
0643 include 'epos.incsem'
0644
0645 conbmxdif=0.
0646 b1=0.
0647 bmax=7.
0648 iomegasave=iomega
0649 iomega=2
0650
0651 eps=1.e-5
0652 pmax=1.d0-dexp(-om1intbc(b1))
0653 pdiff=facdif
0654 if(pmax.gt.eps)then
0655 conbmxdif=drootom(pdiff,pmax,bmax,eps)
0656 endif
0657 iomega=iomegasave
0658
0659 return
0660
0661 end
0662
0663
0664 subroutine conre
0665
0666
0667
0668 include "epos.incems"
0669 include 'epos.inc'
0670
0671 call utpri('conre ',ish,ishini,6)
0672
0673
0674
0675 la=laproj
0676 ma=iabs(maproj)
0677 las=0
0678 mas=0
0679 do l=1,ma
0680 if(la.lt.0)then
0681 if(iabs(idproj).lt.20)then
0682 id=idproj
0683 else
0684 ia=iabs(idproj/10)
0685 is=idproj/iabs(idproj)
0686 if(ia.ne.111.and.ia.ne.222.and.ia.ne.333)id=idproj/10*10
0687 if(ia.eq.111.or. ia.eq.222.or. ia.eq.333)id=idproj/10*10+is
0688 if(ia.eq.213)id=1230*is
0689 endif
0690 else
0691 id=1220
0692 if(rangen().le.(la-las)*1./(ma-mas))id=1120
0693 if(id.eq.1120)las=las+1
0694 mas=mas+1
0695 endif
0696 ic1=idtrai(1,id,1)
0697 ic2=idtrai(2,id,1)
0698 icproj(1,l)=ic1
0699 icproj(2,l)=ic2
0700 enddo
0701
0702
0703
0704 la=latarg
0705 ma=iabs(matarg)
0706 las=0
0707 mas=0
0708 do l=1,ma
0709 if(la.lt.0)then
0710 if(iabs(idtarg).lt.20)then
0711 id=idtarg
0712 else
0713 ia=iabs(idtarg/10)
0714 is=idtarg/iabs(idtarg)
0715 if(ia.ne.111.and.ia.ne.222.and.ia.ne.333)id=idtarg/10*10
0716 if(ia.eq.111.or. ia.eq.222.or. ia.eq.333)id=idtarg/10*10+is
0717 if(ia.eq.213)id=1230*is
0718 endif
0719 else
0720 id=1220
0721 if(rangen().le.(la-las)*1./(ma-mas))id=1120
0722 if(id.eq.1120)las=las+1
0723 mas=mas+1
0724 endif
0725 ic1=idtrai(1,id,1)
0726 ic2=idtrai(2,id,1)
0727 ictarg(1,l)=ic1
0728 ictarg(2,l)=ic2
0729 enddo
0730
0731 call utprix('conre ',ish,ishini,6)
0732 return
0733 end
0734
0735
0736 subroutine conrl
0737
0738
0739
0740 include "epos.incems"
0741 common/nucl1/laproj,maproj,latarg,matarg,core,fctrmx
0742 common/hadr2/iomodl,idproj,idtarg,wexcit
0743
0744
0745
0746 la=latarg
0747 ma=iabs(matarg)
0748 las=0
0749 mas=0
0750 do l=1,ma
0751 if(la.lt.0)then
0752 id=idtarg
0753 else
0754 id=1220
0755 if(rangen().le.(la-las)*1./(ma-mas))id=1120
0756 if(id.eq.1120)las=las+1
0757 mas=mas+1
0758 endif
0759 ic1=idtrai(1,id,1)
0760 ic2=idtrai(2,id,1)
0761 ictarg(1,l)=ic1
0762 ictarg(2,l)=ic2
0763 enddo
0764
0765 return
0766 end
0767
0768
0769 subroutine conwr
0770
0771
0772
0773 include "epos.inc"
0774 include "epos.incems"
0775 double precision XA(210,3),XB(210,3),BQGS,BMAXQGS,BMAXNEX,BMINNEX
0776 COMMON /Q_QGSNEX1/ XA,XB,BQGS,BMAXQGS,BMAXNEX,BMINNEX
0777 common/nucl3/phi,bimp
0778 common /cncl/xproj(mamx),yproj(mamx),zproj(mamx)
0779 *,xtarg(mamx),ytarg(mamx),ztarg(mamx)
0780 parameter(iapmax=208)
0781 double precision bqgs2,bmaxqgs2,bmaxnex2,bminnex2,xan,xbn
0782 common /qgsIInex1/xan(iapmax,3),xbn(iapmax,3)
0783 *,bqgs2,bmaxqgs2,bmaxnex2,bminnex2
0784 common/photrans/phoele(4),ebeam
0785 integer ic(2)
0786
0787 call utpri('conwr ',ish,ishini,6)
0788
0789 bx=cos(phi)*bimp
0790 by=sin(phi)*bimp
0791
0792
0793
0794 nptl=0
0795
0796 if(iokoll.eq.1)then ! precisely matarg collisions
0797
0798 nptl=nptl+1
0799 do 3 i=1,4
0800 3 xorptl(i,nptl)=0
0801 tivptl(1,nptl)=-ainfin
0802 tivptl(2,nptl)=0
0803 istptl(nptl)=1
0804 iorptl(nptl)=-1
0805 jorptl(nptl)=0
0806 do 1 k=1,koll
0807 nptl=nptl+1
0808 do 4 i=1,4
0809 4 xorptl(i,nptl)=0
0810 tivptl(1,nptl)=-ainfin
0811 tivptl(2,nptl)=0
0812 istptl(nptl)=1
0813 iorptl(nptl)=-1
0814 jorptl(nptl)=0
0815 1 continue
0816
0817 elseif(iappl.ne.7)then
0818
0819
0820
0821 do 6 i=1,maproj
0822 nptl=nptl+1
0823 istptl(nptl)=0
0824 iorptl(nptl)=0
0825 jorptl(nptl)=0
0826 if(model.eq.2)then !QGSJet
0827 xproj(i)=XA(i,1)
0828 yproj(i)=XA(i,2)
0829 zproj(i)=XA(i,3)
0830 istptl(nptl)=1
0831 iorptl(nptl)=-1
0832 elseif(model.eq.7.or.model.eq.11)then !QGSJetII
0833 xproj(i)=xan(i,1)
0834 yproj(i)=xan(i,2)
0835 zproj(i)=xan(i,3)
0836 istptl(nptl)=1
0837 iorptl(nptl)=-1
0838 elseif(model.ge.3)then !Gheisha, ...
0839 istptl(nptl)=1
0840 iorptl(nptl)=-1
0841 endif
0842 xorptl(1,nptl)=xproj(i)+bx/2
0843 xorptl(2,nptl)=yproj(i)+by/2
0844 xorptl(3,nptl)=zproj(i)
0845 xorptl(4,nptl)=0
0846 tivptl(1,nptl)=-ainfin
0847
0848
0849 tivptl(2,nptl)= ainfin
0850
0851 6 continue
0852
0853 do 7 i=1,matarg
0854 nptl=nptl+1
0855 istptl(nptl)=0
0856 iorptl(nptl)=0
0857 jorptl(nptl)=0
0858 if(model.eq.2)then !QGSJet
0859 xtarg(i)=XB(i,1)
0860 ytarg(i)=XB(i,2)
0861 ztarg(i)=XB(i,3)
0862 istptl(nptl)=1
0863 iorptl(nptl)=-1
0864 elseif(model.eq.7.or.model.eq.11)then !QGSJetII
0865 xtarg(i)=xbn(i,1)
0866 ytarg(i)=xbn(i,2)
0867 ztarg(i)=xbn(i,3)
0868 istptl(nptl)=1
0869 iorptl(nptl)=-1
0870 elseif(model.ge.3)then !Gheisha, ...
0871 istptl(nptl)=1
0872 iorptl(nptl)=-1
0873 endif
0874 xorptl(1,nptl)=xtarg(i)-bx/2
0875 xorptl(2,nptl)=ytarg(i)-by/2
0876 xorptl(3,nptl)=ztarg(i)
0877 xorptl(4,nptl)=0
0878 tivptl(1,nptl)=-ainfin
0879
0880
0881 tivptl(2,nptl)= ainfin
0882
0883 7 continue
0884 if(abs(idprojin).eq.12)then !electron for fake DIS
0885
0886 nptl=nptl+1
0887 istptl(nptl)=41
0888 iorptl(nptl)=-1
0889 jorptl(nptl)=-1
0890 iorptl(1)=nptl !pi0 (porjectile) coming from lepton
0891 xorptl(1,nptl)=bx/2
0892 xorptl(2,nptl)=by/2
0893 xorptl(3,nptl)=0.
0894 xorptl(4,nptl)=0.
0895 tivptl(1,nptl)=-ainfin
0896 tivptl(2,nptl)=0.
0897
0898 do i=1,matarg
0899 nptl=nptl+1
0900 istptl(nptl)=41
0901 iorptl(nptl)=-1
0902 jorptl(nptl)=-1
0903 xorptl(1,nptl)=xtarg(i)-bx/2
0904 xorptl(2,nptl)=ytarg(i)-by/2
0905 xorptl(3,nptl)=ztarg(i)
0906 xorptl(4,nptl)=0
0907 tivptl(1,nptl)=-ainfin
0908
0909
0910 tivptl(2,nptl)= ainfin
0911 enddo
0912
0913 nptl=nptl+1
0914 istptl(nptl)=0
0915 iorptl(nptl)=maproj+matarg+1
0916 jorptl(nptl)=-1
0917 xorptl(1,nptl)=bx/2
0918 xorptl(2,nptl)=by/2
0919 xorptl(3,nptl)=0.
0920 xorptl(4,nptl)=0.
0921 tivptl(1,nptl)=0.
0922 tivptl(2,nptl)= ainfin
0923 endif
0924
0925 endif
0926
0927 nptl=0
0928 if(iappl.le.2)then
0929 do i=1,maproj
0930 nptl=nptl+1
0931 ic(1)=icproj(1,i)
0932 ic(2)=icproj(2,i)
0933 id=idtra(ic,0,0,0)
0934
0935 call idmass(id,ams)
0936 idptl(nptl)=id
0937 pptl(1,nptl)=0.
0938 pptl(2,nptl)=0.
0939 pptl(3,nptl)=pnullx
0940 pptl(4,nptl)=sqrt(pnullx**2+ams**2)
0941 pptl(5,nptl)=ams
0942 ifrptl(1,nptl)=0
0943 ifrptl(2,nptl)=0
0944 ityptl(nptl)=0
0945 enddo
0946 endif
0947 if(iappl.ne.7)then
0948 do i=1,matarg
0949 nptl=nptl+1
0950 ic(1)=ictarg(1,i)
0951 ic(2)=ictarg(2,i)
0952 id=idtra(ic,0,0,0)
0953
0954 call idmass(id,ams)
0955 idptl(nptl)=id
0956 pptl(1,nptl)=0.
0957 pptl(2,nptl)=0.
0958 pptl(3,nptl)=-pnullx
0959 pptl(4,nptl)=sqrt(pnullx**2+ams**2)
0960 pptl(5,nptl)=ams
0961 ifrptl(1,nptl)=0
0962 ifrptl(2,nptl)=0
0963 ityptl(nptl)=0
0964 enddo
0965 if(abs(idprojin).eq.12)then !electron for fake DIS
0966
0967 nptl=nptl+1
0968 id=idprojin
0969 call idmass(id,ams)
0970 idptl(nptl)=id
0971 pptl(1,nptl)=0.
0972 pptl(2,nptl)=0.
0973 pptl(3,nptl)=sqrt(max(0.,(elepti+ams)*(elepti-ams)))
0974 pptl(4,nptl)=elepti
0975 pptl(5,nptl)=ams
0976 ifrptl(1,nptl)=1
0977 ifrptl(2,nptl)=1
0978 ityptl(nptl)=40
0979
0980 do i=1,matarg
0981 nptl=nptl+1
0982 idptl(nptl)=idptl(maproj+i)
0983 pptl(1,nptl)=0.
0984 pptl(2,nptl)=0.
0985 pptl(3,nptl)=-pnll
0986 pptl(4,nptl)=ebeam
0987 pptl(5,nptl)=pptl(5,maproj+i)
0988 ifrptl(1,nptl)=maproj+i
0989 ifrptl(2,nptl)=maproj+i
0990 ityptl(nptl)=50
0991 enddo
0992
0993 nptl=nptl+1
0994 idptl(nptl)=id
0995 pptl(1,nptl)=phoele(1)
0996 pptl(2,nptl)=phoele(2)
0997 pptl(3,nptl)=phoele(3)
0998 pptl(4,nptl)=phoele(4)
0999 pptl(5,nptl)=ams
1000 ifrptl(1,nptl)=0
1001 ifrptl(2,nptl)=0
1002 ityptl(nptl)=40
1003 endif
1004
1005 else
1006
1007 nptl=nptl+1
1008 id=idproj
1009 call idmass(id,ams)
1010 idptl(nptl)=id
1011 pptl(1,nptl)=0.
1012 pptl(2,nptl)=0.
1013 pptl(3,nptl)=pnullx
1014 pptl(4,nptl)=sqrt(pnullx**2+ams**2)
1015 pptl(5,nptl)=ams
1016 ifrptl(1,nptl)=0
1017 ifrptl(2,nptl)=0
1018 ityptl(nptl)=0
1019 iorptl(nptl)=-1
1020 jorptl(nptl)=0
1021 istptl(nptl)=0
1022 do 5 i=1,4
1023 5 xorptl(i,nptl)=0
1024 tivptl(1,nptl)=0
1025 tivptl(2,nptl)=0
1026 endif
1027
1028
1029
1030
1031 call utprix('conwr ',ish,ishini,6)
1032 return
1033 end
1034
1035
1036 subroutine conxyz(ch,n,x,y,z,ynuc)
1037
1038 include 'epos.inc'
1039
1040 real x(n),y(n),z(n)
1041 character ch*1
1042
1043 massnr=0
1044 iii=0
1045 if(ch.eq.'p')then
1046 massnr=maproj
1047 iii=1
1048 elseif(ch.eq.'t')then
1049 massnr=matarg
1050 iii=2
1051 else
1052 call utstop('conxyz: nucleus neither proj nor targ&')
1053 endif
1054
1055 if(massnr.eq.0)return
1056 if(massnr.gt.n)call utstop('conxyz: massnr.gt.n&')
1057 if(massnr.eq.1)then
1058 x(1)=0
1059 y(1)=0
1060 z(1)=0
1061 return
1062 endif
1063
1064 rad=radnuc(massnr)
1065
1066 if(massnr.ge.10)then !---wood-saxon density---
1067
1068 rad=rad/difnuc(massnr)
1069 cr1=1.+3./rad+6./rad**2+6./rad**3
1070 cr2=3./rad
1071 cr3=3./rad+6./rad**2
1072 do i=1,massnr
1073 1 zuk=rangen()*cr1-1.
1074 if(zuk.le.0.)then
1075 tt=rad*(rangen()**.3333-1.)
1076 elseif(zuk.le.cr2 )then
1077 tt=-log(rangen())
1078 elseif(zuk.lt.cr3 )then
1079 tt=-log(rangen())-log(rangen())
1080 else
1081 tt=-log(rangen())-log(rangen())-log(rangen())
1082 endif
1083 if(rangen().gt.1./(1.+exp(-abs(tt))))goto 1
1084 rim=tt+rad
1085 zz=rim*(2.*rangen()-1.)
1086 rim=sqrt(rim*rim-zz*zz)
1087 z(i)=zz*difnuc(massnr)
1088 call pscs(c,s)
1089 x(i)=rim*c*difnuc(massnr)
1090 y(i)=rim*s*difnuc(massnr)
1091 enddo
1092
1093 elseif(massnr.ge.3)then ! ---gaussian density---
1094
1095 rad=rad*sqrt(2.*massnr/(massnr-1.)) !van hove simulation
1096 do l=1,3
1097 summ=0.
1098 do i=1,massnr-1
1099 j=massnr-i
1100 aks=rad *(rangen()+rangen()+rangen()-1.5)
1101 k=j+1
1102 if(l.eq.1)x(k)=summ-aks*sqrt(float(j)/k)
1103 if(l.eq.2)y(k)=summ-aks*sqrt(float(j)/k)
1104 if(l.eq.3)z(k)=summ-aks*sqrt(float(j)/k)
1105 summ=summ+aks/sqrt(float(j*k))
1106 enddo
1107 if(l.eq.1)x(1)=summ
1108 if(l.eq.2)y(1)=summ
1109 if(l.eq.3)z(1)=summ
1110 enddo
1111
1112 elseif(massnr.eq.2)then ! ---deuteron---
1113
1114 !.........r=t*difnuc(massnr), t~exp(-2*t)*(1-exp(-a*t))
1115 a=radnuc(massnr)
1116 2 t=-0.5*alog(rangen()) !~exp(-2*t)
1117 if(rangen().gt.(1-exp(-a*t))**2)goto2
1118 r=t*difnuc(massnr)
1119 zz=r*(2.*rangen()-1.)
1120 call pscs(c,s)
1121 rxy=sqrt(r*r-zz*zz)
1122 z(1)=0.5*zz
1123 x(1)=0.5*rxy*c
1124 y(1)=0.5*rxy*s
1125 z(2)=-z(1)
1126 x(2)=-x(1)
1127 y(2)=-y(1)
1128
1129 else
1130
1131 stop'conxyz: wrong massnr. '
1132
1133 endif
1134
1135
1136
1137 rmax=(radnuc(massnr)+3)
1138 drnucl(iii)=rmax/mxnucl
1139 nrnucl(iii)=nrnucl(iii)+1
1140 do i=1,massnr
1141 r=sqrt(x(i)**2+y(i)**2+z(i)**2)
1142 k=1+int(r/drnucl(iii))
1143 if(k.le.mxnucl)rnucl(k,iii)=rnucl(k,iii)+1
1144 enddo
1145
1146
1147
1148 do i=1,massnr
1149 z(i)=z(i)/cosh(ynuc)
1150 enddo
1151
1152 return
1153 end
1154
1155
1156 subroutine conini
1157
1158 include 'epos.inc'
1159
1160 imax=max(maproj,matarg)
1161 if(idtargin.eq.0)imax=max(imax,40)
1162 do massnr=1,mamxx
1163 dif=0.54
1164 rad=0.
1165 if(massnr.gt.imax.or.massnr.eq.1)then
1166 dif=0
1167 rad=0
1168 elseif(massnr.eq.197)then
1169 dif=0.562
1170 rad=6.5
1171 elseif(massnr.ge.10)then
1172 rad=1.12*massnr**0.33333-0.86*massnr**(-0.33333)
1173 elseif(massnr.ge.3)then
1174 rad=.9*float(massnr)**.3333
1175 elseif(massnr.eq.2)then
1176 dif=4.316
1177 rad=4.68
1178 endif
1179 difnuc(massnr)=dif
1180 radnuc(massnr)=rad
1181 enddo
1182
1183 end
1184
1185
1186 subroutine xConNuclDens(iii)
1187
1188
1189
1190
1191 include 'epos.inc'
1192 if(model.ne.1)return
1193 massnr=1
1194 if(iii.eq.1)then
1195 massnr=maproj
1196 elseif(iii.eq.2)then
1197 massnr=matarg
1198 endif
1199 if(massnr.eq.1)return
1200 a=1./4.316
1201 b=4.68
1202 write(ifhi,'(a)') '!-----------------------------------------'
1203 write(ifhi,'(a)') '! nuclear density '
1204 write(ifhi,'(a)') '!-----------------------------------------'
1205 write(ifhi,'(a)') 'openhisto'
1206 if(massnr.ge.10)write(ifhi,'(a)')'htyp lin xmod lin ymod lin'
1207 if(massnr.lt.10)write(ifhi,'(a)')'htyp lin xmod lin ymod log'
1208 write(ifhi,'(a)') 'text 0 0 "title nuclear density"'
1209 write(ifhi,'(a)') 'text 0.99 -0.15 " r (fm) "'
1210 write(ifhi,'(a)') 'text 0 0 "yaxis rho(r)"'
1211 write(ifhi,'(a,2e11.3)')'xrange',0.,mxnucl*drnucl(iii)
1212 write(ifhi,'(3a)')'yrange',' 0 ',' auto'
1213 write(ifhi,'(a)') 'array 2'
1214 do j=1,mxnucl
1215 r=(j-0.5)*drnucl(iii)
1216 d=0.5*drnucl(iii)
1217 write(ifhi,'(2e12.4)') r,rnucl(j,iii)/nrnucl(iii)/
1218 * (4./3.*pi*((r+d)**3-(r-d)**3))
1219 enddo
1220 write(ifhi,'(a)') ' endarray'
1221 write(ifhi,'(a)') 'closehisto plot 0-'
1222 write(ifhi,'(a)') 'openhisto'
1223 write(ifhi,'(a)') 'htyp lbo '
1224 write(ifhi,'(a)') 'array 2'
1225 do j=1,mxnucl
1226 r=(j-0.5)*drnucl(iii)
1227 rr=2*r
1228 rho=1
1229 if(massnr.eq.2)then
1230 rho=1.00*((1-exp(-b*a*rr))*exp(-a*rr)/rr)**2
1231 elseif(massnr.eq.197)then
1232 rho=0.16/(1+exp((r-6.5)/0.562))
1233 elseif(massnr.ge.10)then
1234 rho=0.16/(1+exp((r-radnuc(massnr))/difnuc(massnr)))
1235 endif
1236 write(ifhi,'(2e12.4)') r,rho
1237 enddo
1238 write(ifhi,'(a)') ' endarray'
1239 write(ifhi,'(a)') 'closehisto plot 0'
1240 end
1241
1242
1243 subroutine xConThick(iii)
1244
1245 ! plots sigma_pp *T_A (b) (=average number of collisions)
1246 ! T_A = thickness function
1247 ! iii = 1 (proj) or 2 (targ)
1248 !----------------------------------------------------------------
1249 include 'epos.inc'
1250 parameter(iconimax=20,iconkmax=100)
1251 real thick(2,0:iconimax)
1252 imax=iconimax
1253 kmax=iconkmax
1254 if(model.ne.1)return
1255 ramx=mxnucl*drnucl(iii)
1256 do i=0,imax
1257 x=i/float(imax)*ramx
1258 sum=0
1259 rho0=conrho(iii,0.)
1260 h=ramx/kmax
1261 do k=1,kmax
1262 z=k*h
1263 r=sqrt(x*x+z*z)
1264 rho2=conrho(iii,r)
1265 z=(k-0.5)*h
1266 r=sqrt(x*x+z*z)
1267 rho1=conrho(iii,r)
1268 sum=sum+h/6.*(rho0+4*rho1+rho2)
1269 rho0=rho2
1270 enddo
1271 sum=sum*2 ! integral fro -infty to +infty
1272 thick(1,i)=x
1273 thick(2,i)=sum
1274 enddo
1275 write(ifhi,'(a)') 'openhisto'
1276 write(ifhi,'(a)') 'htyp lbo '
1277 write(ifhi,'(a)') 'txt "xaxis b (fm)" '
1278 write(ifhi,'(a)') 'txt "yaxis [s]?pp! T?A! (b) " '
1279 write(ifhi,'(a)') 'array 2'
1280 do i=0,imax
1281 write(ifhi,'(2e12.4)') thick(1,i),sigine/10.*thick(2,i)
1282 enddo
1283 write(ifhi,'(a)') ' endarray'
1284 write(ifhi,'(a)') 'closehisto plot 0'
1285
1286 end
1287
1288
1289 function conrho(iii,r)
1290
1291 ! nuclear density
1292 ! iii = 1 (proj) or 2 (targ)
1293 !----------------------------------------------------------------
1294 include 'epos.inc'
1295 conrho=1.
1296 if(model.ne.1)return
1297 massnr=1
1298 if(iii.eq.1)then
1299 massnr=maproj
1300 elseif(iii.eq.2)then
1301 massnr=matarg
1302 endif
1303 if(massnr.eq.1)return
1304 a=1./4.316
1305 b=4.68
1306 rr=2*r
1307 rho=1
1308 if(massnr.eq.2.and.rr.gt.0.)then
1309 rho=1.00*((1-exp(-b*a*rr))*exp(-a*rr)/rr)**2
1310 elseif(massnr.eq.197)then
1311 rho=0.16/(1+exp((r-6.5)/0.562))
1312 elseif(massnr.ge.10)then
1313 rho=0.16/(1+exp((r-radnuc(massnr))/difnuc(massnr)))
1314 endif
1315 conrho=rho
1316 end
1317
1318
1319
1320
1321
1322
1323
