c Version 2.0, Nov. 1, 1994. program gaussmaker implicit none double precision mass1,mass2,gammav,gamma,rapid,temp double precision pt1,pt2,rx,ry,rz,tau,lambda,b,bmin,bmax,delb double precision bw0,bcut,x,y,z,t,e,px,py,pz,znew,gauss integer*4 IDUM,IDUM2,IY,IV(32) common/RAN2JUNK/IDUM,IDUM2,IY,IV double precision ran2,randy integer ib,nb,i,ident,npts1,npts2,iseed c This program generates phase space points for use by 'cmain.f' c It generates gaussian distributions in time and space and thermal c distributions for the momenta. The phase space points are then c boosted along the z-axis according to the rapidity. The user enters c the temperature, Rx,Ry,Rz and lifetime tau of the source. All c variables correspond to the conv. g(x) = e{-x^2/(2*R**2)}... c The program will run much more efficiently for heavy particles, c (T << mass) if the gaussian routines are used rather than c the subroutine 'relpfind' for generating momenta which are relativistic c in the frame of the source. c The program will generate a fraction of the points, (1-dsqrt(lambda)), c such that correlationss would give a corresponding intercept of c 1+lambda when doing pure Bose problems. c c I think it is easy to figure out the program without comments, c but then, I wrote them. - Scott c --------- Usually you need to edit these lines appropriately ------ c Make ident 0 if particles are not identical, 1 if they are identical ident=0 c Choose the masses of the two particles c mass1=938.28d0 c mass1=939.6d0 c mass2=939.6d0 mass1=139.58d0 c mass2=493.646d0 c mass1=1875.585d0 c mass2=3727.323d0 mass2=mass1 c ----------------------------------------------------- IDUM=-1234 c Choose the number of impact parameters. This is chosen different c than unity only for the purpose of testing the program. open(20,form='formatted',file='phasedescription.dat', + status='unknown') open(21,form='formatted',file='phasespace1.dat',status='unknown') if(ident.eq.0) open(22,form='formatted',file='phasespace2.dat', + status='unknown') nb=1 write(20,*) nb write(6,*) 'what is the rapidity of the thermal source?' read(5,*) rapid gammav=sinh(rapid) gamma=cosh(rapid) write(6,*) 'what is the temperature?' read(5,*) temp pt1=dsqrt(mass1*temp) pt2=dsqrt(mass2*temp) write(6,*) 'what are rx,ry,rz and tau?' read(5,*) rx,ry,rz,tau write(6,*) 'what is lambda?' read(5,*) lambda lambda=dsqrt(lambda) delb=1.0 do 111 ib=1,nb b=delb*(.5d0+dfloat(ib-1)) bmin=b-.5d0*delb bmax=b+.5d0*delb if(nb.eq.1) b=0.d0 npts1=10000 npts2=10000 bw0=1.d0 bcut=1.d0 write(20,*) b,bmin,bmax,npts1,npts2,bw0,bcut do 727 i=1,npts1 c The next 3 lines are for non-rel mom. distributions. px=pt1*gauss() py=pt1*gauss() pz=pt1*gauss() c Use the next line for relativistic thermal distributions. c If the temp. is small (<