EPOS/epos167/epos-con-161.f

   1 c-----------------------------------------------------------------------
   2       subroutine conaa(iret)
   3 c-----------------------------------------------------------------------
   4 c  determines interaction configuration
   5 c-----------------------------------------------------------------------
   6
   7       include 'epos.inc'
   8       include 'epos.incems'
   9       include 'epos.incsem'
  10
  11       common/geom/rmproj,rmtarg,bmax,bkmx
  12       common/nucl3/phi,bimp
  13       common /cncl/xproj(mamx),yproj(mamx),zproj(mamx)
  14      *            ,xtarg(mamx),ytarg(mamx),ztarg(mamx)
  15       common/cncl3/iactpr(mamx),iacttg(mamx)
  16       common/cfacmss/facmss
  17
  18       call utpri('conaa ',ish,ishini,4)
  19
  20       iret=0
  21
  22 c     initialisations
  23 c     ---------------
  24
  25       vel=tanh(ypjtl-yhaha)+tanh(yhaha)
  26
  27 c     determine phi, bimp, coll, iproj, itarg, x/y/zproj, x/y/ztarg
  28 c     ---------------------------------------------------------------
  29
  30            if(iokoll.eq.1)then
  31
  32       koll=matarg
  33       do k=1,koll
  34         do n=1,4
  35           coord(n,k)=0.
  36         enddo
  37       enddo
  38       bimp=0
  39       phi=0
  40       do k=1,koll
  41       bij=bkmx*sqrt(rangen())
  42       bk(k)=bij
  43       iproj(k)=1
  44       itarg(k)=k
  45       enddo
  46       xproj(1)=0
  47       yproj(1)=0
  48       zproj(1)=0
  49       iactpr(1)=1
  50       do k=1,koll
  51       phi=2.*pi*rangen()
  52       xtarg(k)=bij*cos(phi)
  53       ytarg(k)=bij*sin(phi)
  54       ztarg(k)=0
  55       iacttg(k)=1
  56       enddo
  57       lproj(1)=koll
  58       lproj3(1)=0
  59       do k=1,koll
  60       ltarg(k)=1
  61       kproj(1,k)=k
  62       ktarg(k,1)=k
  63       ltarg3(k)=0
  64       if(iscreen.ne.0)then
  65         lproj3(1)=lproj3(1)+1
  66         ltarg3(k)=1
  67         kproj3(1,lproj3(1))=k
  68         ktarg3(k,1)=k
  69       endif
  70       enddo
  71
  72            elseif(maproj.eq.1.and.matarg.eq.1)then
  73
  74       b1=bminim
  75       b2=amin1(bkmx,bmaxim)
  76       if(b1.gt.b2)call utstop('conaa: bmin > bmax&')
  77       bimp=sqrt(b1*b1+(b2*b2-b1*b1)*rangen())
  78       koll=1
  79       do n=1,4
  80         coord(n,1)=0.
  81       enddo
  82       bk(1)=bimp
  83       iproj(1)=1
  84       itarg(1)=1
  85       phi=2.*pi*rangen()
  86       xproj(1)=bimp*cos(phi)
  87       yproj(1)=bimp*sin(phi)
  88       zproj(1)=0
  89       xtarg(1)=0
  90       ytarg(1)=0
  91       ztarg(1)=0
  92       lproj(1)=1
  93       ltarg(1)=1
  94       if(iscreen.ne.0)then
  95         lproj3(1)=1
  96         ltarg3(1)=1
  97       else
  98         lproj3(1)=0
  99         ltarg3(1)=0
 100       endif
 101       kproj(1,1)=1
 102       ktarg(1,1)=1
 103       kproj3(1,1)=1
 104       ktarg3(1,1)=1
 105       iactpr(1)=1  !int(rangen()+1./chad(iclpro)) !active proj.
 106       iacttg(1)=1  !int(rangen()+1./chad(icltar)) !active targ.
 107       if(iactpr(1).eq.0.and.iacttg(1).eq.0)iret=-1
 108
 109            else
 110
 111       call conxyz('p',mamx,xproj,yproj,zproj,ypjtl-yhaha)
 112       call conxyz('t',mamx,xtarg,ytarg,ztarg,yhaha)
 113
 114       bx=0
 115       by=0
 116       if(maproj.gt.0)then
 117       if(bimevt.lt.0)then
 118         b1=bminim
 119         b2=amin1(rmproj+rmtarg,bmaxim)
 120         if(b1.gt.b2)call utstop('conaa: bmin > bmax&')
 121         bimp=sqrt(b1**2+(b2**2-b1**2)*rangen())
 122         if(nbarray.gt.0)bimp=barray(mod(nrevt,nbarray)+1)
 123         if(jpsi.gt.0)then
 124           bimp=b1+(b2-b1)*(float(mod(nrevt,12))+rangen())/12.
 125           bimevt=bimp
 126         endif
 127         phi=phimin+rangen()*(phimax-phimin)
 128       else
 129         phi=phievt
 130         bimp=bimevt
 131       endif
 132       bx=cos(phi)*bimp
 133       by=sin(phi)*bimp
 134       endif
 135       if(jpsi.lt.0)then !modify b
 136         bx=xtarg(1)
 137         by=ytarg(1)
 138       endif
 139       if(maproj.eq.0)goto1000
 140 c      wactv=1./chad(2)              !probability to be "active"
 141       koll=0
 142       do i=1,maproj
 143       lproj(i)=0
 144       lproj3(i)=0
 145       iactpr(i)=1 !int(rangen()+wactv) !active proj. nucleons
 146       enddo
 147       do j=1,matarg
 148       ltarg(j)=0
 149       ltarg3(j)=0
 150       iacttg(j)=1 !int(rangen()+wactv) !active targ. nucleons
 151       enddo
 152       do 12 i=1,maproj
 153       if(iactpr(i).eq.0)goto 12
 154       do 11 j=1,matarg
 155       if(iacttg(j).eq.0)goto 11
 156       if(jpsi.lt.0.and.ztarg(j).le.ztarg(1))goto11
 157       bij=sqrt((xproj(i)+bx-xtarg(j))**2+(yproj(i)+by-ytarg(j))**2)
 158       if(ish.ge.7)write(ifch,*)'i_p:',i,' i_t:',j,' b_ij:',bij
 159       if(bij.gt.bkmx)goto 11
 160
 161       koll=koll+1
 162       if(koll.gt.kollmx)call utstop('conaa: kollmx too small&')
 163       bk(koll)=bij
 164       bkx(koll)=xproj(i)+bx-xtarg(j)
 165       bky(koll)=yproj(i)+by-ytarg(j)
 166       iproj(koll)=i
 167       itarg(koll)=j
 168       lproj(i)=lproj(i)+1
 169       ltarg(j)=ltarg(j)+1
 170       kproj(i,lproj(i))=koll
 171       ktarg(j,ltarg(j))=koll
 172       if(iscreen.ne.0)then
 173         if(bij.gt.bkmxndif)goto 11
 174         lproj3(i)=lproj3(i)+1
 175         ltarg3(j)=ltarg3(j)+1
 176         kproj3(i,lproj3(i))=koll
 177         ktarg3(j,ltarg3(j))=koll
 178       endif
 179
 180 11    continue
 181 12    continue
 182
 183       do k=1,koll
 184         do n=1,4
 185           coord(n,k)=0.
 186         enddo
 187       enddo
 188
 189       if(((iactpr(1).eq.0.and.maproj.eq.1)
 190      &.or.(iacttg(1).eq.0.and.matarg.eq.1))
 191      &.and.iret.eq.0)iret=-1
 192
 193           endif
 194
 195       if(ish.ge.3)write(ifch,*)'koll=',koll
 196       if(koll.eq.0)goto 1001
 197
 198
 199 c     determine coord
 200 c     ---------------
 201       do kl=1,koll
 202       i=iproj(kl)
 203       j=itarg(kl)
 204       dist=ztarg(j)-zproj(i)
 205       coord(1,kl)=(xproj(i)+xtarg(j))*0.5
 206       coord(2,kl)=(yproj(i)+ytarg(j))*0.5
 207       coord(3,kl)=(zproj(i)+ztarg(j))*0.5
 208       coord(4,kl)=dist/vel
 209       enddo
 210
 211
 212
 213 c     exit
 214 c     ----
 215 1000  continue
 216       if(ish.ge.5)then
 217       write(ifch,*)'ia,x/y/zproj:'
 218       do mm=1,maproj
 219       write(ifch,*)mm,iactpr(mm),xproj(mm),yproj(mm),zproj(mm)
 220       enddo
 221       write(ifch,*)'ia,x/y/ztarg:'
 222       do mm=1,matarg
 223       write(ifch,*)mm,iacttg(mm),xtarg(mm),ytarg(mm),ztarg(mm)
 224       enddo
 225       write(ifch,*)'iret',iret
 226       endif
 227       call utprix('conaa ',ish,ishini,4)
 228       return
 229
 230 1001  continue !iret=1 causes redo of whole collision
 231       iret=1
 232       if(ish.ge.3)then
 233       write(ifch,*)
 234       write(ifch,*)'***** subroutine conaa:'
 235       write(ifch,*)'***** no nucleon pair found --> no interaction'
 236       write(ifch,*)
 237       endif
 238       goto 1000
 239
 240       end
 241
 242 c-----------------------------------------------------------------------
 243       subroutine xGeometry(iii)
 244 c-----------------------------------------------------------------------
 245       include 'epos.inc'
 246       include 'epos.incems'
 247       include 'epos.incsem'
 248       common/xgeom/nnn,naa(kollmx),nbb(kollmx)
 249       character*5 fmt1,fmt2
 250
 251       if(iii.eq.0)then
 252
 253       do k=1,kollmx
 254         naa(k)=0
 255       enddo
 256       nnn=0
 257
 258
 259       elseif(iii.eq.1)then
 260
 261       ngl=0
 262       do k=1,koll
 263         r=bk(k)
 264         if(r.le.sqrt(sigine/10./pi))ngl=ngl+1
 265       enddo
 266       if(ngl.ne.0)then
 267         nnn=nnn+1
 268         naa(ngl)=naa(ngl)+1
 269       endif
 270
 271       elseif(iii.eq.2)then
 272
 273       if(xpar1.eq.0..and.xpar2.eq.0.)then
 274        print*,'---------- xGeometry -----------------------'
 275        return
 276       endif
 277       x1=1-0.01*xpar2
 278       x2=1-0.01*xpar1
 279       kmx=0
 280       nbb(1)=naa(1)
 281       do k=2,kollmx
 282         if(naa(k).ne.0.)kmx=k
 283         nbb(k)=nbb(k-1)+naa(k)
 284       enddo
 285       k1=0
 286       k2=0
 287       do k=1,kmx
 288         x=nbb(k)/float(nnn)
 289         if(x.lt.x1)k1=k
 290         if(x.lt.x2)k2=k
 291       enddo
 292       k1=k1+1
 293       k2=k2+1
 294       x=0
 295       av=0
 296       su=0
 297       do k=1,kmx
 298         xb=x
 299         x=nbb(k)/float(nnn)
 300         y=naa(k)/float(nnn)
 301         dx=x-xb
 302         p=0.
 303         if(k.eq.k1)then
 304           p=(x-x1)/dx
 305           p1=p
 306         elseif(k.eq.k2)then
 307           p=(x2-xb)/dx
 308           p2=p
 309         elseif(k.gt.k1.and.k.lt.k2)then
 310           p=1
 311         endif
 312       av=av+y*p*k
 313       su=su+y*p
 314       enddo
 315       av=av/su
 316       n1=nint(100*p1)
 317       n2=nint(100*p2)
 318       if(n1.eq.0)then
 319         k1=k1+1
 320         n1=100
 321       endif
 322       if(k1.le.9)fmt1='i1,4x'
 323       if(k1.gt.9.and.k1.le.99)fmt1='i2,3x'
 324       if(k1.gt.99.and.k1.le.999)fmt1='i3,2x'
 325       if(k1.gt.999.and.k1.le.9999)fmt1='i4,1x'
 326       if(k2.le.9)fmt2='i1,4x'
 327       if(k2.gt.9.and.k2.le.99)fmt2='i2,3x'
 328       if(k2.gt.99.and.k2.le.999)fmt2='i3,2x'
 329       if(k2.gt.999.and.k2.le.9999)fmt2='i4,1x'
 330       write(6,'(i4,a,i5,a,'//fmt1//',i6,a,i5,a,'//fmt2//',5x,a,f8.2)')
 331      &       nint(xpar2),':MIN',n1,'%',k1
 332      &      ,nint(xpar1),':MAX',n2,'%',k2   ,'av:',av
 333       endif
 334
 335       end
 336
 337 c-----------------------------------------------------------------------
 338       function conbmx()
 339 c-----------------------------------------------------------------------
 340       double precision om1intbc,p,eps
 341       include 'epos.inc'
 342       include 'epos.incsem'
 343
 344       conbmx=0.
 345       b1=0.
 346       b2=7.
 347       eps=5.0d-3
 348       p=1.d0-dexp(-om1intbc(b2))
 349       if(p.gt.2.d0*eps)return
 350
 351       ntry=0
 352
 353 10    ntry=ntry+1
 354       b=b1+.5*(b2-b1)
 355
 356       p=(1.d0-dexp(-om1intbc(b)))
 357
 358       if(p.gt.eps)then
 359        if(p.gt.2.d0*eps)then
 360         b1=b
 361        else
 362         conbmx=b
 363        return
 364        endif
 365       else
 366        if(p.lt.eps/5.d0)then
 367         b2=b
 368        else
 369         conbmx=b
 370         return
 371        endif
 372       endif
 373
 374       if(ntry.le.1000)goto 10
 375       write(ifmt,*)'Too much try in conbmx ... bmax=',b
 376       conbmx=b
 377       return
 378
 379       end
 380
 381 c-----------------------------------------------------------------------
 382       function conbmxndif()
 383 c-----------------------------------------------------------------------
 384       double precision om1intbc,p,eps
 385       include 'epos.inc'
 386       include 'epos.incsem'
 387
 388
 389       iomegasave=iomega
 390       iomega=2
 391       conbmxndif=0.
 392       b1=0.
 393       b2=7.
 394       bkmxndif=b2
 395       eps=zbcut
 396       p=1.d0-dexp(-om1intbc(b2))
 397       if(p.gt.2.d0*eps)goto 100
 398
 399       ntry=0
 400
 401 10    ntry=ntry+1
 402       b=b1+.5*(b2-b1)
 403
 404       p=(1.d0-dexp(-om1intbc(b)))
 405
 406       if(p.gt.eps)then
 407        if(p.gt.2.d0*eps)then
 408         b1=b
 409        else
 410         conbmxndif=b
 411        goto 100
 412        endif
 413       else
 414        if(p.lt.eps/5.d0)then
 415         b2=b
 416        else
 417         conbmxndif=b
 418         goto 100
 419        endif
 420       endif
 421
 422       if(ntry.le.1000)goto 10
 423       write(ifmt,*)'Too much try in conbmxndif ... bkmxndif=',b
 424       conbmxndif=b
 425  100  iomega=iomegasave
 426       return
 427
 428       end
 429
 430 c-----------------------------------------------------------------------
 431       function conbmxdif()
 432 c-----------------------------------------------------------------------
 433 c find b to have (1-exp(-om))pmax=pdiff
 434 c-----------------------------------------------------------------------
 435       double precision om1intbc,pmax,drootom,pdiff
 436       include 'epos.inc'
 437       include 'epos.incsem'
 438
 439       conbmxdif=0.
 440       b1=0.
 441       bmax=7.
 442       iomegasave=iomega
 443       iomega=2
 444
 445       eps=1.e-5
 446       pmax=1.d0-dexp(-om1intbc(b1))
 447       pdiff=facdif
 448       if(pmax.gt.eps)then
 449         conbmxdif=drootom(pdiff,pmax,bmax,eps)
 450       endif
 451       iomega=iomegasave
 452
 453       return
 454
 455       end
 456
 457 c-----------------------------------------------------------------------
 458       subroutine conre
 459 c-----------------------------------------------------------------------
 460 c  initializes remnants
 461 c-----------------------------------------------------------------------
 462       include "epos.incems"
 463       include 'epos.inc'
 464
 465       call utpri('conre ',ish,ishini,6)
 466
 467 c     proj
 468 c     ----
 469       la=laproj
 470       ma=iabs(maproj)
 471       las=0
 472       mas=0
 473       do l=1,ma
 474       if(la.lt.0)then
 475       ia=iabs(idproj/10)
 476       is=idproj/iabs(idproj)
 477       if(ia.ne.111.and.ia.ne.222.and.ia.ne.333)id=idproj/10*10
 478       if(ia.eq.111.or. ia.eq.222.or. ia.eq.333)id=idproj/10*10+is
 479       if(ia.eq.213)id=1230*is
 480       if(iabs(idproj).lt.20)id=idproj
 481       else
 482       id=1220
 483       if(rangen().le.(la-las)*1./(ma-mas))id=1120
 484       if(id.eq.1120)las=las+1
 485       mas=mas+1
 486       endif
 487       ic1=idtrai(1,id,1)
 488       ic2=idtrai(2,id,1)
 489       icproj(1,l)=ic1
 490       icproj(2,l)=ic2
 491       enddo
 492
 493 c     targ
 494 c     ----
 495       la=latarg
 496       ma=iabs(matarg)
 497       las=0
 498       mas=0
 499       do l=1,ma
 500       if(la.lt.0)then
 501       ia=iabs(idtarg/10)
 502       is=idtarg/iabs(idtarg)
 503       if(ia.ne.111.and.ia.ne.222.and.ia.ne.333)id=idtarg/10*10
 504       if(ia.eq.111.or. ia.eq.222.or. ia.eq.333)id=idtarg/10*10+is
 505       if(ia.eq.213)id=1230*is
 506       if(iabs(idtarg).lt.20)id=idtarg
 507       else
 508       id=1220
 509       if(rangen().le.(la-las)*1./(ma-mas))id=1120
 510       if(id.eq.1120)las=las+1
 511       mas=mas+1
 512       endif
 513       ic1=idtrai(1,id,1)
 514       ic2=idtrai(2,id,1)
 515       ictarg(1,l)=ic1
 516       ictarg(2,l)=ic2
 517       enddo
 518
 519       call utprix('conre ',ish,ishini,6)
 520       return
 521       end
 522
 523 c-----------------------------------------------------------------------
 524       subroutine conrl
 525 c-----------------------------------------------------------------------
 526 c  initializes target remnant in case of appl lepton
 527 c-----------------------------------------------------------------------
 528       include "epos.incems"
 529       common/nucl1/laproj,maproj,latarg,matarg,core,fctrmx
 530       common/hadr2/iomodl,idproj,idtarg,wexcit
 531
 532 c     targ
 533 c     ----
 534       la=latarg
 535       ma=iabs(matarg)
 536       las=0
 537       mas=0
 538       do l=1,ma
 539       if(la.lt.0)then
 540       id=idtarg
 541       else
 542       id=1220
 543       if(rangen().le.(la-las)*1./(ma-mas))id=1120
 544       if(id.eq.1120)las=las+1
 545       mas=mas+1
 546       endif
 547       ic1=idtrai(1,id,1)
 548       ic2=idtrai(2,id,1)
 549       ictarg(1,l)=ic1
 550       ictarg(2,l)=ic2
 551       enddo
 552
 553       return
 554       end
 555
 556 c-----------------------------------------------------------------------
 557       subroutine conwr
 558 c-----------------------------------------------------------------------
 559 c     writes /cptl/
 560 c-----------------------------------------------------------------------
 561       include "epos.inc"
 562       include "epos.incems"
 563       double precision XA(64,3),XB(64,3),BQGS,BMAXQGS,BMAXNEX,BMINNEX
 564       COMMON /Q_QGSNEX1/ XA,XB,BQGS,BMAXQGS,BMAXNEX,BMINNEX
 565       common/nucl3/phi,bimp
 566       common /cncl/xproj(mamx),yproj(mamx),zproj(mamx)
 567      *,xtarg(mamx),ytarg(mamx),ztarg(mamx)
 568       parameter(iapmax=207)
 569       double precision bqgs2,bmaxqgs2,bmaxnex2,bminnex2,xan,xbn
 570       common /qgsIInex1/xan(iapmax,3),xbn(iapmax,3)
 571      *,bqgs2,bmaxqgs2,bmaxnex2,bminnex2
 572       integer ic(2)
 573
 574       call utpri('conwr ',ish,ishini,6)
 575
 576       bx=cos(phi)*bimp
 577       by=sin(phi)*bimp
 578
 579 c     write /cptl/
 580 c     ------------
 581       nptl=0
 582
 583            if(iokoll.eq.1)then   ! precisely matarg collisions
 584
 585       nptl=nptl+1
 586       do 3 i=1,4
 587 3     xorptl(i,nptl)=0
 588       tivptl(1,nptl)=-ainfin
 589       tivptl(2,nptl)=0
 590       istptl(nptl)=1
 591       iorptl(nptl)=-1
 592       jorptl(nptl)=0
 593       do 1 k=1,koll
 594       nptl=nptl+1
 595       do 4 i=1,4
 596 4     xorptl(i,nptl)=0
 597       tivptl(1,nptl)=-ainfin
 598       tivptl(2,nptl)=0
 599       istptl(nptl)=1
 600       iorptl(nptl)=-1
 601       jorptl(nptl)=0
 602 1     continue
 603
 604            elseif(iappl.ne.7)then
 605
 606       do 6 i=1,maproj
 607       nptl=nptl+1
 608       istptl(nptl)=0
 609       iorptl(nptl)=0
 610       jorptl(nptl)=0
 611       if(model.eq.2)then       !QGSJet
 612       xproj(i)=XA(i,1)
 613       yproj(i)=XA(i,2)
 614       zproj(i)=XA(i,3)
 615       istptl(nptl)=1
 616       iorptl(nptl)=-1
 617       elseif(model.eq.7)then       !QGSJetII
 618       xproj(i)=xan(i,1)
 619       yproj(i)=xan(i,2)
 620       zproj(i)=xan(i,3)
 621       istptl(nptl)=1
 622       iorptl(nptl)=-1
 623       elseif(model.ge.3)then       !Gheisha, ...
 624       istptl(nptl)=1
 625       iorptl(nptl)=-1
 626       endif
 627       xorptl(1,nptl)=xproj(i)+bx/2
 628       xorptl(2,nptl)=yproj(i)+by/2
 629       xorptl(3,nptl)=zproj(i)
 630       xorptl(4,nptl)=0
 631       tivptl(1,nptl)=-ainfin
 632 c for visualisation uncomment
 633 c-c   tivptl(1,nptl)=-100
 634       tivptl(2,nptl)= ainfin
 635 6     continue
 636       do 7 i=1,matarg
 637       nptl=nptl+1
 638       istptl(nptl)=0
 639       iorptl(nptl)=0
 640       jorptl(nptl)=0
 641       if(model.eq.2)then       !QGSJet
 642       xtarg(i)=XB(i,1)
 643       ytarg(i)=XB(i,2)
 644       ztarg(i)=XB(i,3)
 645       istptl(nptl)=1
 646       iorptl(nptl)=-1
 647       elseif(model.eq.7)then       !QGSJetII
 648       xtarg(i)=xbn(i,1)
 649       ytarg(i)=xbn(i,2)
 650       ztarg(i)=xbn(i,3)
 651       istptl(nptl)=1
 652       iorptl(nptl)=-1
 653       elseif(model.ge.3)then       !Gheisha, ...
 654       istptl(nptl)=1
 655       iorptl(nptl)=-1
 656       endif
 657       xorptl(1,nptl)=xtarg(i)-bx/2
 658       xorptl(2,nptl)=ytarg(i)-by/2
 659       xorptl(3,nptl)=ztarg(i)
 660       xorptl(4,nptl)=0
 661       tivptl(1,nptl)=-ainfin
 662 c for visualisation uncomment
 663 c-c   tivptl(1,nptl)=-100
 664       tivptl(2,nptl)= ainfin
 665 7     continue
 666
 667         endif
 668
 669       nptl=0
 670       if(iappl.le.2)then
 671       do i=1,maproj
 672       nptl=nptl+1
 673       ic(1)=icproj(1,i)
 674       ic(2)=icproj(2,i)
 675       id=idtra(ic,0,0,3)
 676       call idmass(id,ams)
 677       idptl(nptl)=id
 678       pptl(1,nptl)=0.
 679       pptl(2,nptl)=0.
 680       pptl(3,nptl)=pnullx
 681       pptl(4,nptl)=sqrt(pnullx**2+ams**2)
 682       pptl(5,nptl)=ams
 683       ifrptl(1,nptl)=0
 684       ifrptl(2,nptl)=0
 685       ityptl(nptl)=0
 686       enddo
 687       endif
 688       if(iappl.ne.7)then
 689       do i=1,matarg
 690       nptl=nptl+1
 691       ic(1)=ictarg(1,i)
 692       ic(2)=ictarg(2,i)
 693       id=idtra(ic,0,0,3)
 694       call idmass(id,ams)
 695       idptl(nptl)=id
 696       pptl(1,nptl)=0.
 697       pptl(2,nptl)=0.
 698       pptl(3,nptl)=-pnullx
 699       pptl(4,nptl)=sqrt(pnullx**2+ams**2)
 700       pptl(5,nptl)=ams
 701       ifrptl(1,nptl)=0
 702       ifrptl(2,nptl)=0
 703       ityptl(nptl)=0
 704       enddo
 705
 706       else
 707
 708       nptl=nptl+1
 709       id=idproj
 710       call idmass(id,ams)
 711       idptl(nptl)=id
 712       pptl(1,nptl)=0.
 713       pptl(2,nptl)=0.
 714       pptl(3,nptl)=pnullx
 715       pptl(4,nptl)=sqrt(pnullx**2+ams**2)
 716       pptl(5,nptl)=ams
 717       ifrptl(1,nptl)=0
 718       ifrptl(2,nptl)=0
 719       ityptl(nptl)=0
 720       iorptl(nptl)=-1
 721       jorptl(nptl)=0
 722       istptl(nptl)=0
 723       do 5 i=1,4
 724  5      xorptl(i,nptl)=0
 725       tivptl(1,nptl)=0
 726       tivptl(2,nptl)=0
 727       endif
 728
 729 c     exit
 730 c     ----
 731
 732       call utprix('conwr ',ish,ishini,6)
 733       return
 734       end
 735
 736 c------------------------------------------------------------------------
 737       subroutine conxyz(ch,n,x,y,z,ynuc)
 738 c-----------------------------------------------------------------------
 739       include 'epos.inc'
 740
 741       real x(n),y(n),z(n)
 742       character ch*1
 743
 744       if(ch.eq.'p')then
 745       massnr=maproj
 746       iii=1
 747       elseif(ch.eq.'t')then
 748       massnr=matarg
 749       iii=2
 750       else
 751       call utstop('conxyz: nucleus neither proj nor targ&')
 752       endif
 753
 754       if(massnr.eq.0)return
 755       if(massnr.gt.n)call utstop('conxyz: massnr.gt.n&')
 756       if(massnr.eq.1)then
 757       x(1)=0
 758       y(1)=0
 759       z(1)=0
 760       return
 761       endif
 762
 763       rad=radnuc(massnr)
 764
 765       if(massnr.ge.10)then !---wood-saxon density---
 766
 767         rad=rad/difnuc(massnr)
 768         cr1=1.+3./rad+6./rad**2+6./rad**3
 769         cr2=3./rad
 770         cr3=3./rad+6./rad**2
 771         do i=1,massnr
 772    1      zuk=rangen()*cr1-1.
 773           if(zuk.le.0.)then
 774             tt=rad*(rangen()**.3333-1.)
 775           elseif(zuk.le.cr2 )then
 776             tt=-log(rangen())
 777           elseif(zuk.lt.cr3 )then
 778             tt=-log(rangen())-log(rangen())
 779           else
 780             tt=-log(rangen())-log(rangen())-log(rangen())
 781           endif
 782           if(rangen().gt.1./(1.+exp(-abs(tt))))goto 1
 783           rim=tt+rad
 784           zz=rim*(2.*rangen()-1.)
 785           rim=sqrt(rim*rim-zz*zz)
 786           z(i)=zz*difnuc(massnr)
 787           call pscs(c,s)
 788           x(i)=rim*c*difnuc(massnr)
 789           y(i)=rim*s*difnuc(massnr)
 790         enddo
 791
 792       elseif(massnr.ge.3)then  ! ---gaussian density---
 793
 794         rad=rad*sqrt(2.*massnr/(massnr-1.))   !van hove simulation
 795         do l=1,3
 796           summ=0.
 797           do i=1,massnr-1
 798             j=massnr-i
 799             aks=rad *(rangen()+rangen()+rangen()-1.5)
 800             k=j+1
 801             if(l.eq.1)x(k)=summ-aks*sqrt(float(j)/k)
 802             if(l.eq.2)y(k)=summ-aks*sqrt(float(j)/k)
 803             if(l.eq.3)z(k)=summ-aks*sqrt(float(j)/k)
 804             summ=summ+aks/sqrt(float(j*k))
 805           enddo
 806           if(l.eq.1)x(1)=summ
 807           if(l.eq.2)y(1)=summ
 808           if(l.eq.3)z(1)=summ
 809         enddo
 810
 811       elseif(massnr.eq.2)then  ! ---deuteron---
 812
 813         !.........r=t*difnuc(massnr), t~exp(-2*t)*(1-exp(-a*t))
 814         a=radnuc(massnr)
 815   2     t=-0.5*alog(rangen())  !~exp(-2*t)
 816         if(rangen().gt.(1-exp(-a*t))**2)goto2
 817         r=t*difnuc(massnr)
 818         zz=r*(2.*rangen()-1.)
 819         call pscs(c,s)
 820         rxy=sqrt(r*r-zz*zz)
 821         z(1)=0.5*zz
 822         x(1)=0.5*rxy*c
 823         y(1)=0.5*rxy*s
 824         z(2)=-z(1)
 825         x(2)=-x(1)
 826         y(2)=-y(1)
 827
 828       else
 829
 830         stop'conxyz: wrong massnr.     '
 831
 832       endif
 833
 834 c...plot preparation
 835
 836       rmax=(radnuc(massnr)+3)
 837       drnucl(iii)=rmax/mxnucl
 838       nrnucl(iii)=nrnucl(iii)+1
 839       do i=1,massnr
 840         r=sqrt(x(i)**2+y(i)**2+z(i)**2)
 841         k=1+int(r/drnucl(iii))
 842         if(k.le.mxnucl)rnucl(k,iii)=rnucl(k,iii)+1
 843       enddo
 844
 845 c...lorentz trafo
 846
 847       do i=1,massnr
 848       z(i)=z(i)/cosh(ynuc)
 849       enddo
 850
 851       return
 852       end
 853
 854 c-----------------------------------------------------------------------
 855       subroutine conini
 856 c-----------------------------------------------------------------------
 857       include 'epos.inc'
 858
 859       imax=max(maproj,matarg)
 860       if(idtargin.eq.0)imax=max(imax,40)
 861       do massnr=1,mamxx
 862         dif=0.54
 863         if(massnr.gt.imax.or.massnr.eq.1)then
 864           dif=0
 865           rad=0
 866         elseif(massnr.eq.197)then
 867           dif=0.562
 868           rad=6.5
 869         elseif(massnr.ge.10)then
 870           rad=1.12*massnr**0.33333-0.86*massnr**(-0.33333)
 871         elseif(massnr.ge.3)then
 872           rad=.9*float(massnr)**.3333
 873         elseif(massnr.eq.2)then
 874           dif=4.316
 875           rad=4.68
 876         endif
 877         difnuc(massnr)=dif
 878         radnuc(massnr)=rad
 879       enddo
 880       end
 881
 882 c-----------------------------------------------------------------------
 883       subroutine xConNuclDens(iii)
 884 c-----------------------------------------------------------------------
 885 c plots distribution of nucleons in nuclei
 886 c  iii = 1 (proj) or 2 (targ)
 887 c-----------------------------------------------------------------------
 888       include 'epos.inc'
 889       if(model.ne.1)return
 890       if(iii.eq.1)then
 891         massnr=maproj
 892       elseif(iii.eq.2)then
 893         massnr=matarg
 894       endif
 895       a=1./4.316
 896       b=4.68
 897       write(ifhi,'(a)') '!-----------------------------------------'
 898       write(ifhi,'(a)') '!          nuclear density          '
 899       write(ifhi,'(a)') '!-----------------------------------------'
 900       write(ifhi,'(a)')       'openhisto'
 901       if(massnr.ge.10)write(ifhi,'(a)')'htyp lin xmod lin ymod lin'
 902       if(massnr.lt.10)write(ifhi,'(a)')'htyp lin xmod lin ymod log'
 903       write(ifhi,'(a)')       'text 0 0 "title nuclear density"'
 904       write(ifhi,'(a)')       'text 0.99 -0.15 " r (fm) "'
 905       write(ifhi,'(a)')       'text 0 0 "yaxis rho(r)"'
 906       write(ifhi,'(a,2e11.3)')'xrange',0.,mxnucl*drnucl(iii)
 907       write(ifhi,'(3a)')'yrange',' 0 ',' auto'
 908       write(ifhi,'(a)')       'array 2'
 909       do j=1,mxnucl
 910       r=(j-0.5)*drnucl(iii)
 911       d=0.5*drnucl(iii)
 912       write(ifhi,'(2e12.4)')  r,rnucl(j,iii)/nrnucl(iii)/
 913      *                     (4./3.*pi*((r+d)**3-(r-d)**3))
 914       enddo
 915       write(ifhi,'(a)')       '  endarray'
 916       write(ifhi,'(a)')       'closehisto plot 0-'
 917       write(ifhi,'(a)')       'openhisto'
 918       write(ifhi,'(a)')       'htyp lbo '
 919       write(ifhi,'(a)')       'array 2'
 920       do j=1,mxnucl
 921       r=(j-0.5)*drnucl(iii)
 922       rr=2*r
 923       rho=1
 924       if(massnr.eq.2)then
 925         rho=1.00*((1-exp(-b*a*rr))*exp(-a*rr)/rr)**2
 926       elseif(massnr.eq.197)then
 927         rho=0.16/(1+exp((r-6.5)/0.562))
 928       elseif(massnr.ge.10)then
 929         rho=0.16/(1+exp((r-radnuc(massnr))/difnuc(massnr)))
 930       endif
 931       write(ifhi,'(2e12.4)')  r,rho
 932       enddo
 933       write(ifhi,'(a)')       '  endarray'
 934       write(ifhi,'(a)')       'closehisto plot 0'
 935       end
 936
 937 c-----------------------------------------------------------------------
 938       subroutine xConThick(iii)
 939 c-----------------------------------------------------------------------
 940       ! plots sigma_pp *T_A (b)  (=average number of collisions)
 941       ! T_A = thickness function
 942       !  iii = 1 (proj) or 2 (targ)
 943       !----------------------------------------------------------------
 944       include 'epos.inc'
 945       parameter(iconimax=20,iconkmax=100)
 946       real thick(2,0:iconimax)
 947       imax=iconimax
 948       kmax=iconkmax
 949       if(model.ne.1)return
 950       ramx=mxnucl*drnucl(iii)
 951       do i=0,imax
 952         x=i/float(imax)*ramx
 953         sum=0
 954         rho0=conrho(iii,0.)
 955         h=ramx/kmax
 956         do k=1,kmax
 957           z=k*h
 958           r=sqrt(x*x+z*z)
 959           rho2=conrho(iii,r)
 960           z=(k-0.5)*h
 961           r=sqrt(x*x+z*z)
 962           rho1=conrho(iii,r)
 963           sum=sum+h/6.*(rho0+4*rho1+rho2)
 964           rho0=rho2
 965         enddo
 966         sum=sum*2 ! integral fro -infty to +infty
 967         thick(1,i)=x
 968         thick(2,i)=sum
 969       enddo
 970       write(ifhi,'(a)')       'openhisto'
 971       write(ifhi,'(a)')       'htyp lbo '
 972       write(ifhi,'(a)')       'txt "xaxis b (fm)" '
 973       write(ifhi,'(a)')       'txt "yaxis [s]?pp! T?A! (b) " '
 974       write(ifhi,'(a)')       'array 2'
 975       do i=0,imax
 976         write(ifhi,'(2e12.4)') thick(1,i),sigine/10.*thick(2,i)
 977       enddo
 978       write(ifhi,'(a)')       '  endarray'
 979       write(ifhi,'(a)')       'closehisto plot 0'
 980
 981       end
 982 c-----------------------------------------------------------------------
 983       function conrho(iii,r)
 984 c-----------------------------------------------------------------------
 985       ! nuclear density
 986       !  iii = 1 (proj) or 2 (targ)
 987       !----------------------------------------------------------------
 988       include 'epos.inc'
 989       conrho=1.
 990       if(model.ne.1)return
 991       if(iii.eq.1)then
 992         massnr=maproj
 993       elseif(iii.eq.2)then
 994         massnr=matarg
 995       endif
 996       a=1./4.316
 997       b=4.68
 998       rr=2*r
 999       rho=1
1000       if(massnr.eq.2)then
1001         rho=1.00*((1-exp(-b*a*rr))*exp(-a*rr)/rr)**2
1002       elseif(massnr.eq.197)then
1003         rho=0.16/(1+exp((r-6.5)/0.562))
1004       elseif(massnr.ge.10)then
1005         rho=0.16/(1+exp((r-radnuc(massnr))/difnuc(massnr)))
1006       endif
1007       conrho=rho
1008       end
1009
1010
1011
1012
1013
1014
1015