EPOS/epos167/epos-yyy-168.f

   1 c----------------------------------------------------------------------
   2       subroutine hnbaaa156(ip,iret)
   3 c----------------------------------------------------------------------
   4 c  microcanonical decay of cluster ip via loop over hnbmet
   5 c----------------------------------------------------------------------
   6       include 'epos.inc'
   7       common/cxyzt/xptl(mxptl),yptl(mxptl),zptl(mxptl),tptl(mxptl)
   8      *,optl(mxptl),uptl(mxptl),sptl(mxptl),rptl(mxptl,3)
   9       parameter(maxp=500)
  10       common/confg/np,amass(maxp),ident(maxp),pcm(5,maxp),wtxlog,wtlog
  11       common/citer/iter,itermx
  12       double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
  13       common /cttaus/  tpro,zpro,ttar,ztar,ttaus,detap,detat
  14       integer jc(nflav,2)
  15       double precision p(5),c(5)
  16       parameter(maxit=50000)
  17       common/count/nacc,nrej,naccit(maxit),nptot,npit(maxit)
  18       dimension be(4),pe(5),pa(5)
  19       common/yradx/yrad(maxp),phirad(maxp)
  20       common/xxxspecsy/ndrop(-4:4,-4:4,-4:4)
  21       common/cdelzet/delzet,delsgr /cvocell/vocell
  22       call utpri('hnbaaa',ish,ishini,4)
  23
  24       if(ish.ge.3)then
  25       write(ifch,140)sngl(ttaus)
  26   140 format(/' ----------------------------------'/
  27      *'    droplet decay  at tau =',f6.2/
  28      *' ----------------------------------')
  29       write(ifch,*)'droplet:'
  30       call alist('&',ip,ip)
  31       endif
  32
  33       iret=0
  34       do j=1,5
  35       c(j)=pptl(j,ip)
  36       enddo
  37
  38       call idquac(ip,nqi,nsi,nai,jc)
  39       keu=jc(1,1)-jc(1,2)
  40       ked=jc(2,1)-jc(2,2)
  41       kes=jc(3,1)-jc(3,2)
  42       kec=jc(4,1)-jc(4,2)
  43       keb=jc(5,1)-jc(5,2)
  44       ket=jc(6,1)-jc(6,2)
  45       !print*,'droplet uds=',keu,ked,kes,'   E=',pptl(5,ip)
  46
  47       volu=4./3.*pi*radptl(ip)**3
  48
  49       if(volu.le.0)call utstop('hnbaaa: volume = 0&')
  50       if(volu.lt.0.01)then
  51         call utmsg('hnbaaa')
  52         write(ifch,*)'*****  very small volume:',volu
  53         write(ifch,*)
  54      *  'id:',idptl(ip),' r:',radptl(ip),' m:',pptl(5,ip)
  55         call utmsgf
  56       endif
  57
  58     !~~~~~redefine energy in case of radial flow
  59       amin=utamnu(keu,ked,kes,kec,keb,ket,4)   !utamnu(...,4) and not utamnu(...,5)
  60       aumin=amuseg                             !otherwise droplet from remnant decay
  61       ipo=ip                                   !could be too light after flow
  62       if(ityptl(ip).eq.60)ipo=iorptl(ip)
  63       tecmor=pptl(5,ipo)
  64       if(iappl.eq.4.or.iorsdf.ne.3
  65      &.or.ityptl(ip).eq.40.or.ityptl(ip).eq.50)then !not for droplets from remnants
  66         yrmax=0
  67       else
  68         yrmax=max(0.0,yradmx+yradmi*alog10(engy/200.))
  69                 if(maproj.eq.1.and.matarg.eq.1)then
  70           yrmax=max(0.0,yradpp+yradpi*alog10(engy/1800.))
  71                   !aumin=amin
  72                   !if(yrmax.gt.0.2)print*,'===',tecmor,aamin,yrmax
  73                 endif
  74       endif
  75       fradflo=1.
  76       if(yrmax.ne.0.)
  77      &fradflo=1./((sinh(yrmax)*yrmax-cosh(yrmax)+1.)/(yrmax**2/2.))
  78       tecm=pptl(5,ip)
  79       tecmxx=tecm
  80       if(yrmax.gt.0..and.tecmor.gt.aumin
  81      &  .and.tecm*fradflo.gt.amin) then
  82         ! redefine energy to account for collective flow
  83         ! \int_0^yrmax f(y) d2y = E_new (effective mass)
  84         ! \int_0^yrmax cosh(y) f(y) d2y = E_old
  85         tecm=tecm*fradflo
  86         if(tecm.lt.amin)stop'aaahnb: small mass. should not happen.   '
  87       else
  88         yrmax=0.
  89       endif
  90
  91     !~~~~~redefine energy in case of long coll flow
  92       if(iappl.eq.4.or.iorsdf.ne.3
  93      &.or.ityptl(ip).eq.40.or.ityptl(ip).eq.50)then !not for droplets from remnants
  94         yco=0
  95       else
  96        if(ylongmx.lt.0.)then
  97         yco=delzet * 1.75
  98        else
  99         yco=ylongmx
 100        endif
 101       endif
 102       tecmx=tecm
 103       if(yco.gt.0..and.tecmor.gt.aumin) then
 104         tecm=tecm/sinh(yco)*yco
 105       else
 106         yco=0.
 107       endif
 108       !print*,'========= cluster energy: ',pptl(5,ip),tecmx,tecm
 109
 110     !~~~~~~~~~redefine volume~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 111
 112 ccc      vocri=tecm/epscri(ioclude)   !??????????????????????????????
 113 ccc      volu=max(vocri,vocell)       !????????????????????????????
 114
 115     !~~~~~~~~~decay~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 116       call hnbini(iret)
 117       !if(iret.ne.0)write(ifch,*)'***** unsucessfull hnbini *****'
 118       if(iret.ne.0)goto1000
 119       if(ioinct.ge.1)goto1
 120
 121       do iter=1,itermx
 122         naccit(iter)=0
 123         call hnbmet
 124       enddo
 125
 126 1     continue
 127
 128       if(ioceau.eq.1.and.iappl.eq.1)call xhnbte(ip)
 129
 130     !~~~~~~~~~~long coll flow -> particles~~~~~~~~~~~~~~~~
 131       if(yco.gt.0.) then
 132         errlim=0.0001
 133         tecm=tecmx
 134  611    energ=0.
 135         do i=1,np
 136           yrad(i)=(2*rangen()-1)*yco
 137           be(3)=sinh(yrad(i))
 138           be(4)=cosh(yrad(i))
 139           energ=energ+be(4)*pcm(4,i)-be(3)*pcm(3,i)
 140         enddo
 141         if(abs(energ-tecm).gt.0.1) goto 611
 142                 !print*,'===== energy after flow boosts',energ,'   soll: ',tecm
 143         do j=1,4
 144                   pe(j)=0.
 145                 enddo
 146         do i=1,np
 147           be(1)= 0
 148           be(2)= 0
 149           be(3)= sinh(yrad(i))
 150           be(4)= cosh(yrad(i))
 151           call utlob3(1,be(1),be(2),be(3),be(4),1e0
 152      *         , pcm(1,i), pcm(2,i), pcm(3,i), pcm(4,i))
 153           do j=1,4
 154           pe(j)=pe(j)+pcm(j,i)
 155           enddo
 156         enddo
 157                 pe(5)=sqrt(pe(4)**2-pe(3)**2-pe(2)**2-pe(1)**2)
 158                 !write(6,'(a,5e11.3)')'flow boosts',pe
 159         do j=1,4
 160                   pa(j)=0.
 161                 enddo
 162         do i=1,np
 163           call utlob3(1,pe(1),pe(2),pe(3),pe(4),pe(5)
 164      *         , pcm(1,i), pcm(2,i), pcm(3,i), pcm(4,i))
 165           do j=1,4
 166                     pa(j)=pa(j)+pcm(j,i)
 167                   enddo
 168         enddo
 169                 pa(5)=sqrt(pa(4)**2-pa(3)**2-pa(2)**2-pa(1)**2)
 170                 !write(6,'(a,5e11.3)')' cms boost ',pa
 171         esoll=tecm
 172         scal=1.
 173         do ipass=1,200
 174           sum=0.
 175           do  j=1,np
 176             do k=1,3
 177               pcm(k,j)=scal*pcm(k,j)
 178             enddo
 179             pcm(4,j)=sqrt(pcm(1,j)**2+pcm(2,j)**2+pcm(3,j)**2
 180      *           +amass(j)**2)
 181             sum=sum+pcm(4,j)
 182           enddo
 183           scal=esoll/sum
 184           !write(6,*)'ipass,scal,e,esoll:'
 185           !    $         ,ipass,scal,sum,esoll
 186           if(abs(scal-1.).le.errlim) goto301
 187         enddo
 188  301    continue
 189         do j=1,4
 190                   pa(j)=0.
 191                 enddo
 192         do i=1,np
 193           do j=1,4
 194                     pa(j)=pa(j)+pcm(j,i)
 195                   enddo
 196         enddo
 197         pa(5)=sqrt(pa(4)**2-pa(3)**2-pa(2)**2-pa(1)**2)
 198         !write(6,'(a,5e11.3)')' rescaling ',pa
 199       endif
 200
 201     !~~~~~~~~~~radial flow -> particles~~~~~~~~~~~~~~~~~~
 202       if(yrmax.gt.0.) then
 203         fecc=0
 204                 aa=1
 205                 bb=0
 206                 cc=0
 207                 dd=1
 208         if(ityptl(ip).eq.60)then
 209                   ipo=iorptl(ip)
 210                   xx=uptl(ipo)   ! <x**2>
 211                   yy=optl(ipo)   ! <y**2>
 212                   xy=desptl(ipo) ! <x*y>
 213                   dta=0.5*abs(xx-yy)
 214                   ev1=(xx+yy)/2+sqrt(dta**2+xy**2)
 215                   ev2=(xx+yy)/2-sqrt(dta**2+xy**2)
 216                   if(xy.lt.0..and.dta.ne.0.)then
 217                     theta=0.5*atan(-xy/dta)
 218                   elseif(xy.gt.0..and.dta.ne.0.)then
 219                     theta=-0.5*atan(xy/dta)
 220                   else
 221                     theta=0
 222                   endif
 223                   !eccx=(yy-xx)/(yy+xx)
 224           yy=ev1
 225                   xx=ev2
 226                   ecc=(yy-xx)/(yy+xx)
 227                   !print*,eccx,ecc,theta
 228                   fecc=facecc*ecc
 229                   fecc=min(0.3,fecc)
 230                   phiclu=phievt+theta
 231           aa=cos(phiclu)
 232           bb=sin(phiclu)
 233                   cc=-sin(phiclu)
 234                   dd=cos(phiclu)
 235         endif
 236                 errlim=0.0001
 237         tecm=tecmxx
 238  610    energ=0.
 239         do i=1,np
 240           yrad(i)=sqrt(rangen())
 241           phirad(i)=2.*pi*rangen()
 242           bex=-dsinh(dble(yrad(i)*yrmax))*cos(phirad(i))*(1+fecc)
 243           bey=-dsinh(dble(yrad(i)*yrmax))*sin(phirad(i))*(1-fecc)
 244           be(1)=aa*bex+cc*bey
 245           be(2)=bb*bex+dd*bey
 246           be(3)=-0d0
 247           be(4)=sqrt(1+be(1)**2+be(2)**2)
 248           bp=0d0
 249           do k=1,3
 250             bp=bp+pcm(k,i)*be(k)
 251           enddo
 252           en=be(4)*pcm(4,i)+bp
 253           energ=energ+en
 254         enddo
 255         if(abs(energ-tecm).gt.0.1) goto 610
 256         energ=0.
 257         do i=1,np
 258           bex=dsinh(dble(yrad(i)*yrmax))*cos(phirad(i))*(1+fecc)
 259           bey=dsinh(dble(yrad(i)*yrmax))*sin(phirad(i))*(1-fecc)
 260           be(1)=aa*bex+cc*bey
 261           be(2)=bb*bex+dd*bey
 262           be(3)=0d0
 263           be(4)=sqrt(1+be(1)**2+be(2)**2)
 264           call utlob3(1,be(1),be(2),be(3),be(4),1e0
 265      *         , pcm(1,i), pcm(2,i), pcm(3,i), pcm(4,i))
 266           energ=energ+pcm(4,i)
 267         enddo
 268         esoll=tecm
 269         scal=1.
 270         do ipass=1,200
 271           sum=0.
 272           do  j=1,np
 273             do k=1,3
 274               pcm(k,j)=scal*pcm(k,j)
 275             enddo
 276             pcm(4,j)=sqrt(pcm(1,j)**2+pcm(2,j)**2+pcm(3,j)**2
 277      *           +amass(j)**2)
 278             sum=sum+pcm(4,j)
 279           enddo
 280           scal=esoll/sum
 281           !write(6,*)'ipass,scal,e,esoll:'
 282           ! $         ,ipass,scal,sum,esoll
 283           if(abs(scal-1.).le.errlim) goto300
 284         enddo
 285  300    continue
 286       else
 287         do n=1,np
 288           yrad(n)=0.
 289           phirad(n)=0.
 290         enddo
 291       endif
 292     !~~~~~~~~~~~~~~~
 293
 294       nptlb=nptl
 295       do n=1,np
 296         nptl=nptl+1
 297         if(nptl.gt.mxptl)call utstop('hnbptl: mxptl too small&')
 298         idptl(nptl)=ident(n)
 299         do j=1,4
 300           p(j)=pcm(j,n)
 301         enddo
 302         p(5)=amass(n)
 303         call utlob2(-1,c(1),c(2),c(3),c(4),c(5),p(1),p(2),p(3),p(4),10)
 304         do j=1,5
 305           pptl(j,nptl)=p(j)
 306         enddo
 307         if(tecmor.gt.aumin)then
 308           ityptl(nptl)=60
 309         else
 310           ityptl(nptl)=19
 311         endif
 312         if(ityptl(ip).eq.60)then
 313           if(ityptl(nptl).eq.60)then
 314             ipo=iorptl(ip)
 315             xx=uptl(ipo)        ! <x**2>
 316             yy=optl(ipo)        ! <y**2>
 317             rini=sqrt(5./3.*(xx+yy)) !<r**2>=3/5*R**2 for sphere of radius R
 318             r=1.15*rini*yrad(n) !yrad=y/ymax
 319             tau=2.25/sqrt(yrad(n)**2+0.04)-0.75
 320             z=xorptl(3,ipo)
 321             t=xorptl(4,ipo)
 322            !zeta=0.5*log((t+z)/(t-z))-0.5*delzet+2*0.5*delzet*rangen()
 323             zeta=0.5*log((p(4)+p(3))/(p(4)-p(3)))
 324             z=tau*sinh(zeta)
 325             t=tau*cosh(zeta)
 326             xorptl(1,nptl)=xorptl(1,ipo)+r*cos(phirad(n))
 327             xorptl(2,nptl)=xorptl(2,ipo)+r*sin(phirad(n))
 328             xorptl(3,nptl)=z
 329             xorptl(4,nptl)=t
 330           else
 331             xorptl(1,nptl)=xorptl(1,ip)
 332             xorptl(2,nptl)=xorptl(2,ip)
 333             xorptl(3,nptl)=xorptl(3,ip)
 334             xorptl(4,nptl)=xorptl(4,ip)
 335           endif
 336         endif
 337       enddo
 338
 339       if(ish.ge.3)then
 340         write(ifch,*)'decay products:'
 341         call alist('&',nptlb+1,nptl)
 342         if(ish.ge.5)then
 343           write(ifch,*)'momentum sum:'
 344           do kk=1,5
 345             pptl(kk,nptl+1)=0
 346             do ii=nptlb+1,nptl
 347               pptl(kk,nptl+1)=pptl(kk,nptl+1)+pptl(kk,ii)
 348             enddo
 349             pptl(kk,nptl+2)=c(kk)
 350           enddo
 351           call alist('&',nptl+1,nptl+2)
 352         endif
 353       endif
 354
 355  1000 continue
 356       call utprix('hnbaaa',ish,ishini,4)
 357       return
 358       end