program main c version 2.0, November 1, 1994 c -------------------------------------------- c version 1.0 Dec. 7, 1992 c version 1.1 June 28, 1993 (significant cleaning, several errors c found) c version 1.11 August 9, 1993 (Max Potekhin's bug corrected) c Version 1.2 Dec. 11, 1993 (Minor changes in how coul. arrays are c treated and addition of Breit-Wigner type correlations.) c Version 1.21 Jan. 4, 1994 (Minor corrections to how Breit-Wigner c information is passed, and inclusion of p-alpha and p-pi+. If this c version of cmain is used corresponding version of corrstarter must c be used as well.) c c In version 2, smearing is included along to account for detector c resoloution along with a more transparent explanation of various c direction-defining conventions. The main program and the subroutines c 'transport' have been changed significantly. The routines c 'decompose' and 'smear' are new. See directions for explanations. c In version 2.01 random number generator was changed to ran2 c c ______________________________________________ c NOTES ON UNITS c All momenta and energy are in MeV/c - MeV and all distances and times c are in fermi or fm/c. c ______________________________________________ implicit none c Parameters are set and explained in "cparams.for" include 'cparams.f' c ______________________________________________ c These are the phase space coordinates to be read off c file 21 and 22. c The second index refers to the impact parameter. double precision xx1(nphase,nba),yy1(nphase,nba),zz1(nphase,nba), + tt1(nphase,nba),ppx1(nphase,nba),ppy1(nphase,nba), + ppz1(nphase,nba),ee1(nphase,nba) double precision xx2(nphase,nba),yy2(nphase,nba),zz2(nphase,nba), + tt2(nphase,nba),ppx2(nphase,nba),ppy2(nphase,nba), + ppz2(nphase,nba),ee2(nphase,nba) double precision x1,x2,y1,y2,z1,z2,t1,t2 double precision px1,px2,py1,py2,pz1,pz2,e1,e2 c These are the smeared momenta double precision psx1,psy1,psz1,psx2,psy2,psz2 c ---------------------------------------------- c These are used if the denominator is to be calculated with c iden equal to 2. double precision pppx1(npmax),pppx2(npmax), + pppy1(npmax),pppy2(npmax), + pppz1(npmax),pppz2(npmax),eee1(npmax),eee2(npmax) integer iibb1(npmax),iibb2(npmax),jjj1(npmax),jjj2(npmax) double precision www1(npmax),www2(npmax),w1,w2,wbar1,wbar2 c ______________________________________________ c If the two particles are not identical, get rid of this c equivalence statement equivalence (xx1,xx2),(yy1,yy2),(zz1,zz2),(tt1,tt2), + (ppx1,ppx2),(ppy1,ppy2),(ppz1,ppz2),(ee1,ee2),(np1,np2) c ______________________________________________ c Below are reduced mass(rmass), a constant pdotk that depends c on the mass, symm. weights, and the momentum arrays. c 'msame' is unity if the masses are the same, zero otherwise double precision rmass,wsym1,wsym2,wanti1,wanti2,amom(nka), + akk(nka),mass1,mass2 integer msame c ---------------------------------------------- c First line below, The number of partial waves which are corrected c for strong interaction, number of mom. bins, c ang. momenta of partial waves: 'iq' is z1*z2 and iconvpass c equals the parameter idirconv, which needs to be passed to a sub, c and 'vs' which is the passing of 'v_source'. double precision vs integer npart,nkmax,l1,l2 integer iq,iconvpass c ______________________________________________ c These refer to then number phase space points to be read off c units 21 and 22. They also refer to info. about the c inpact parameters used in weighting. These are decscribed below. integer npts1(nba),npts2(nba),np1(nba),np2(nba),nb double precision bw0(nba),bcut(nba),b(nba),bmin(nba),bmax(nba) double precision bweight(nba),bwtot c ______________________________________________ c 'corr' is the function that calculates the correlation c functions. double precision dcorr,corr c ______________________________________________ c The cN(nka) are the correlation functions. The dN(nka) are the c counters of uncorrelated pairs (the denominator of the corr- c elation functions). The ciaN(nka) are the correlation c due to impact parameters. double precision c0(nka),c1(nka),c2(nka),c3(nka) double precision d0(nka),d1(nka),d2(nka), + d3(nka),cia0(nka),cia1(nka),cia2(nka),cia3(nka) c ______________________________________________ c The c3d(nka) are the correlation functions in 3d space. c The d3d(nka) are the counters of uncorrelated pairs c (the denominator of the corr- c elation functions). The ciaN(nka) are the correlation c due to impact parameters. double precision c3d(nka,nka,nka),d3d(nka,nka,nka) double precision cia3d(nka,nka,nka) double precision delk3d integer ix3d,iy3d,iz3d c ______________________________________________ c These are various projections of the momentum. c 'kx,ky,kz' are the projections of the relative momentum c defined in the routine "decompose". c The relative momentum is defined by (m1*p2-m2*p1)/(m1+m2), c with the zeroth component projected out in the frame of the c total momentum. 'mmom' is the vector sum of these three c components. ptot2 and mom2 are the total and c momentum squared. 'pdotk' is dot product of the total and c relative momentum. c 'mom' is the invariant relative momentum in the frame of the pair. c 'mom2' is -mom**2 and 'mtest' is used to test wheter mom2 c is small enough to store in the correlation arrays. c 'p1dotp2' is p1 dot p2. 'kdotr' is the relative momentum dotted c into the relative spatial coord. c in the frame of the pair. 'r' is the magnitude of the relative c position of the two particles in the frame of the pair. double precision kx,ky,kz,kdotr,mom,mtest,mom2,mmom double precision ptot2,pdotk,p1dotp2,r c ______________________________________________ c 'phi' is an azimuthal angle, 'pi'=3.14..., 'ci' is dsqrt(-1). double precision phi,pi parameter (pi=3.14159265358979323844d0) double complex ci c ______________________________________________ c 'eta' and 'gamow' are the sommerfeld parameter and gamow c factor used in coulomb wave functions. c If you wish to gamow correct the answer simply divide the c correlation functions by gamow at the end of the program. c Otherwise you can leave those references commented out double precision gamow double precision eta c ______________________________________________ c These are counters integer ib,ib1,ib2,iibb integer j1,j2,itry1,itry2 integer ii,kk,iii c ______________________________________________ c 'ia' and 'ifilter' are basically logicals used to signal c acceptance. integer ifilter,ia c ______________________________________________ c 'ibb' is used to test whether the calculation will be c for one impact parameter, b=0, which allows rotating the c points due to cylindrical symmetry. integer ibb c ______________________________________________ c 'nmax' is the number of pairs used to calc. the c correlation function. 1E6 makes a good function. integer nmax c ______________________________________________ c These are used for the random number generator 'ran2' which c comes from Num. recipes integer*4 IDUM,IDUM2,IY,IV(32) common/RAN2JUNK/IDUM,IDUM2,IY,IV double precision randy,ran2 c ______________________________________________ c These are used to describe the relative wave function. double complex cdel1(1000,nka),cdel2(1000,nka), + c1far(nka),c2far(nka) double precision phase1,phase2 double complex c1phase,c2phase,cgamma common /spwave/cdel1,cdel2,c1far,c2far c ______________________________________________ c These are used in the coulomb wave function routines. double complex acouly(nka),acfact1(nka),acfact2(nka), + achype(40,nka),achype1(5,nka),achype2(5,nka) common /couljunk/acouly,acfact1,acfact2,achype,achype1,achype2 c ______________________________________________ c These are used by the function corr. For a description see c below. common /corrinfo/rmass,wsym1,wsym2,wanti1,wanti2,iq, + npart,l1,l2,mass1,mass2 c ---------------------------------------------- c These next variables are used if one wishes to approximate the c wave function squared with a gaussian and Breit-wigner combination. c Most of the variables used by this have "cheap" in the label c to signify this is a cheap way to get an answer. c The necessary things to enter is the energy of the resonance, c (as measured from threshhold) 'echeap', the full width of the c resonance 'gcheap' ( which one can include the experimental c resolution if one wants and the answer will be correct), c the total ang. mom of the resonance 'jcheap', the c spins of the two particles 'j1cheap' and 'j2cheap', c the branching ration into the channel 'bcheap' c the size of the state 'rcheap' which is rather arbitrary, c (But if it is too small statistics will suffer and if it c is too big the answer won't make sense. choose <= one fermi), c and finally a parameter 'icheap' which determines whether one c ignores this part or not. If 'icheap' equals zero one ignores c this stuff, if 'icheap' equals 1, then one goes forward with c this part of the calc. 'iicheap' equals 'icheap', only 'icheap' c is not a parameter so it can be passed to subs. One can do this c part of the calc. and still do the coulomb and partial wave c corrections as well. integer icheap,ijcheap common /cheapint/icheap,ijcheap double precision jcheap,j1cheap,j2cheap,echeap, + gcheap,rcheap,bcheap,momcheap common /cheapreal/jcheap,j1cheap,j2cheap, + echeap,gcheap,rcheap,bcheap,momcheap c --------------------------------------------- c unit 20 has information regarding the impact parameters and c number of phase space points stored. The phase space points c are recorded on units 21 and 22. c unit 25 is initialized by the program corrstarter and carries c information about particle types, masses,etc. along with the c relative wave functions and the relative-momentum-mesh used c to plot the correlation function. open(unit=20,file='phasedescription.dat',status='old', + form='formatted') open(unit=25,status='old',form='unformatted') gamow=1.d0 icheap=iicheap iconvpass=idirconv ci=(0.d0,1.d0) c Initialize random number generator by using negative number for IDUM IDUM=-12345 c write(6,*) 'what is idum? (enter a negative number)' c read(5,*) idum c c Here we read information from unit 25 which was created by the c initialization program. The init. program informs this program c by telling what kind of particles we are working with. c 'rmass' is the reduced mass. c 'wsym1' and 'wsym2' are the weights of symmetric pieces. c If there is only one s-channel then wsym2 will be zero. c 'wanti1' and 'wanti2' are the weights of the anti-symmetric c pieces. c wanti1+wanti2+wsym1+wsym2=1.0 c 'iq' is the product of the charges. c 'npart' is the number of partial waves for which the c Schr. eq. is solved. (incl. strong interactions) c 'l1' and 'l2' are the angular momenta of the partial waves. c cdel1 and cdel2 are corrections to the partial waves c due to strong interactions. c c write(6,*) 'how many points to calculate?' read(5,*) nmax read(25) mass1,mass2 rmass=mass1*mass2/(mass1+mass2) pdotk=mass2**2-mass1**2 msame=0 if(dabs(mass1-mass2).lt.0.0001d0) msame=1 read(25) wsym1,wsym2,wanti1,wanti2 read(25) iq,ijcheap read(25) npart,l1,l2 read(25) nkmax if(nkmax.gt.nka) write(6,*) 'redimension arrays, nka too small' do 990 kk=1,nkmax read(25) amom(kk),akk(kk) if(npart.ge.1) read(25) (cdel1(ii,kk),ii=1,1000) if(npart.ge.2) read(25) (cdel2(ii,kk),ii=1,1000) phase1=0.d0 if(npart.ge.1) read(25) phase1 phase2=0.d0 if(npart.ge.2) read(25) phase2 c c0 - c3 refer to the numerator for the correlation function c cia0 - cia3 are the same thing but with phi**2=1, such that c they will calculate that part of the correlation due to impact c parameter averaging. c d0 -d3 refer to denominators. The "0" suffix refers to c including all angles, while the "1" - "3" suffix refer to c the three different angles. It does not cost you anything c computationally to include the angle cuts. c0(kk)=0.d0 c1(kk)=0.d0 c2(kk)=0.d0 c3(kk)=0.d0 cia0(kk)=0.d0 cia1(kk)=0.d0 cia2(kk)=0.d0 cia3(kk)=0.d0 d0(kk)=0.d0 d1(kk)=0.d0 d2(kk)=0.d0 d3(kk)=0.d0 eta=0.d0 if(iq.ne.0) eta=dfloat(iq)*(rmass/137.036d0)/amom(kk) c -----these are used for assymptotic wave functions---- c1phase=cgamma(1.d0+dfloat(l1)+ci*eta) c2phase=cgamma(1.d0+dfloat(l2)+ci*eta) c1phase=c1phase/dconjg(c1phase) c2phase=c2phase/dconjg(c2phase) c1far(kk)=(cdexp(2.d0*ci*phase1)-1.d0)*c1phase/2.d0 c2far(kk)=(cdexp(2.d0*ci*phase2)-1.d0)*c2phase/2.d0 c ______________________________________________ c Note that we are calculating some numbers used in the c coulomb wave function routines. Here, I have chosen to c calculate the correlation function on a mesh, meaning for c certain values of the momentum. This greatly speeds the c calculation of the coul. wave function(a factor of 10), but could c easily be changed if computers get so fast that the speed c increase is irrelavant. Also, if one wishes to calc.the c correction to the strong int. without using a momentum mesh, c it would really take forever. if(iq.ne.0) call coulset(kk,amom(kk),eta) do 991 ii=1,1000 if(npart.lt.2) cdel2(ii,kk)=(0.d0,0.d0) if(npart.eq.0) cdel1(ii,kk)=(0.d0,0.d0) 991 continue 990 continue if(icheap.eq.1) call cheapset() mtest=-akk(nkmax)**2 c Now get all the phase space pt.s from units 21 and 22 c First read unit 20 to find out how many pt.s to read. c nb is the no. of impact parameters done. c 'b' refers to the impact parameter, 'bmin' and 'bmax' c refer to the range of impact parameters this file is supposed c to represent. This is needed to do the weighting correctly. c 'npts' is the number of phase space pt.s to read in for c that impact parameter c 'bw0 is the weighting factor for that impact parameter. c Choose it such that bw0(ib)*npts(ib) is prop. to the number of c particles emitted for a single collsion for that impact c parameter. c 'np(ib)' is the number of phase space pt.s that c passed the filter. Only these are recorded in the arrays. c bcut(ib) is the probability that the event passes some c impact parameter cut. This cut must be independent of c viewing two particles. If there is no extra cut (aside from c viewing two particles) then each value of bcut should be unity. read(20,*) nb if(nb.gt.nba) write(6,*) + 'too many impact parameters. nba is too small' open(unit=21,form='formatted',status='old', + file='phasespace1.dat') if(ident.ne.1) open(unit=22,form='formatted',status='old', + file='phasespace2.dat') c The next two lines give an example of how to open a file of a c different name. c write(infilename,'(a8,i2.2,a4)') 'dataname',ib,'.dat' c open(unit=21,form='formatted',status='old',file=infilename) do 111 ib=1,nb read(20,*) b(ib),bmin(ib),bmax(ib), + npts1(ib),npts2(ib),bw0(ib),bcut(ib) c Check to see if this calculation is only working with 'b'=0. c if this is true, we can increase statistics by rotating the c points later on in the program, using the spherical symmetry. c Look for the variable 'ibb' later in the program to see the c extra rotation. ibb=0 if(nb.eq.1.and.dabs(b(1)).lt.0.001d0) ibb=1 c *********************************************** c Read in phase space pt.s for first particle np1(ib)=0 do 222 j1=1,npts1(ib) read(21,*,END=222) px1,py1,pz1,x1,y1,z1,t1 c filter1 is the one-particle filter that can not depend on c the azimuthal angle phi, because the points will be rotated c about the beam axis for efficiency. call filter1(mass1,px1,py1,pz1,ifilter) if(ifilter.eq.1) then np1(ib)=np1(ib)+1 if(np1(ib).gt.nphase) write(6,*) 'nphase too small' ppx1(np1(ib),ib)=px1 ppy1(np1(ib),ib)=py1 ppz1(np1(ib),ib)=pz1 ee1(np1(ib),ib)= + dsqrt(px1**2+py1**2+pz1**2+mass1**2) xx1(np1(ib),ib)=x1 yy1(np1(ib),ib)=y1 zz1(np1(ib),ib)=z1 tt1(np1(ib),ib)=t1 endif 222 continue c *********************************************** c Read in phase space pt.s for second particle. c *********************************************** c Skip reading in the second set if the particles are identical. if(ident.eq.1) go to 335 np2(ib)=0 do 333 j1=1,npts2(ib) read(22,*,END=333) px2,py2,pz2,x2,y2,z2,t2 call filter1(mass2,px2,py2,pz2,ifilter) if(ifilter.eq.1) then np2(ib)=np2(ib)+1 if(np2(ib).gt.nphase) write(6,*) 'nphase too small' ppx2(np2(ib),ib)=px2 ppy2(np2(ib),ib)=py2 ppz2(np2(ib),ib)=pz2 ee2(np2(ib),ib)= + dsqrt(px2**2+py2**2+pz2**2+mass2**2) xx2(np2(ib),ib)=x2 yy2(np2(ib),ib)=y2 zz2(np2(ib),ib)=z2 tt2(np2(ib),ib)=t2 endif 333 continue 335 write(6,742) b(ib),mass1,np1(ib),mass2,np2(ib) 742 format(1x,'b,m1,np1,m2,np2=',f8.3,2x,f9.3,2x, + i7,2x,f9.3,2x,i7) 111 continue close(21) if(ident.ne.1) close(22) c *********************************************** write(6,*) 'phase space points read in' c *********************************************** c In the next section of the code we calculate the numerator part c of the correlation function. c ______________________________________________ c First generate the weights for choosing different impact c parameters. The weights are cummulative ones. bwtot=0.d0 do 112 ib=1,nb bweight(ib)=bwtot+dfloat(np1(ib))*dfloat(np2(ib))* + (bw0(ib)**2)*(bmax(ib)**2-bmin(ib)**2)*bcut(ib) bwtot=bweight(ib) 112 continue do 113 ib=1,nb bweight(ib)=bweight(ib)/bwtot 113 continue c ______________________________________________ c itry1=0 iii=0 wbar1=0.d0 wbar2=0.d0 do 727 ii=1,nmax c Pick phase space points, itry1 is used to see how often a pair c works 444 itry1=itry1+1 c First pick the impact parameter randy=ran2() do 449 iibb=1,nb if(randy.lt.bweight(iibb)) then ib=iibb go to 448 endif 449 continue c Now choose the phase space points 448 j1=1+int(dfloat(np1(ib))*ran2()) j2=1+int(dfloat(np2(ib))*ran2()) if(ident.eq.1.and.j2.eq.j1) go to 444 x1=xx1(j1,ib) y1=yy1(j1,ib) z1=zz1(j1,ib) t1=tt1(j1,ib) px1=ppx1(j1,ib) py1=ppy1(j1,ib) pz1=ppz1(j1,ib) e1=ee1(j1,ib) x2=xx2(j2,ib) y2=yy2(j2,ib) z2=zz2(j2,ib) t2=tt2(j2,ib) px2=ppx2(j2,ib) py2=ppy2(j2,ib) if(ibb.eq.1) then phi=2.d0*pi*ran2() call ranrot(px2,py2,x2,y2,phi) endif pz2=ppz2(j2,ib) e2=ee2(j2,ib) c test momentum first. If it is too big go back c and pick another point. These next lines are redundant with what c is done in 'transform' but for wide acceptances, where one does c not wish to calc. pairs with large rel. momenta, it makes the c calc. more efficient to quickly screen out pairs with relative c mometa larger than dsqrt(-mtest). c 'mom' will be the magnitude of the rel. momentum in the frame c where the total momentum is zero. For identical mass part.s c this is the invariant momentum. c p1dotp2=e2*e1-px2*px1-py2*py1-pz2*pz1 mom2=0.25d0*(mass1**2+mass2**2-2.d0*p1dotp2) if(msame.eq.0) then ptot2=mass1**2+mass2**2+2.d0*p1dotp2 mom2=mom2-(0.25d0*pdotk*pdotk/ptot2) endif if(mom2.lt.mtest.and.iden.eq.1) go to 444 if(mom2.lt.mtest.and.iden.eq.2.and.iii.ge.npmax) go to 444 if(mom2.gt.0.d0) write(6,*) 'A: -mom2>0, =',-mom2 if(mom2.gt.0.d0) go to 444 mom=dsqrt(dabs(mom2)) c Now randomly rotate the points about the beam axis. c If 'filter2' does not have any phi dependence,this c upcoming rotation can be skipped. phi=2.d0*pi*ran2() call ranrot(px1,py1,x1,y1,phi) call ranrot(px2,py2,x2,y2,phi) c Now see if these two points are both measured by the apparatus. call filter2(mass1,px1,py1,pz1,e1, + mass2,px2,py2,pz2,e2,ifilter) if(ifilter.eq.0) go to 444 c Now find the relative position r and kdotr, the dot-product of c the relative momentum and the relative position. c call transform(mass1,px1,py1,pz1,e1, + x1,y1,z1,t1,mass2,px2,py2,pz2,e2,x2,y2,z2,t2,mom, + r,kdotr) c Find out if the rel. mom. is less than akk(nkmax), otherwise repick. c Now find the value of kk corresponding to the momentum call kkfind(kk,mom,akk,nkmax) c Now do four angle cuts. If the angle cuts are positive, ia=1, c and the value of corr is added to the correlation function. c c0(kk) - c3(kk) represent the numerator of the correlation c function. cia0(kk) - cia3(kk) represent the correlation due c to impact parameter averaging. The correlation function for c a single set of pt.s is the wave function squared and is found c by the function corr. if(kk.le.nkmax) kdotr=kdotr*amom(kk)/mom dcorr=1.d0 gamow=1.d0 if(r.lt.1.0E-5) go to 444 if(r.lt.1.0E5.and.kk.le.nkmax) + dcorr=corr(kk,amom(kk),kdotr,r) if(kk.le.nkmax.and.igamow.eq.1) then eta=dfloat(iq)*(rmass/137.036d0)/amom(kk) gamow=2*pi*eta/(exp(2*pi*eta)-1.d0) dcorr=dcorr/gamow endif c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c If you will construct your denom. using only two-part. events, c then you need the next several lines of code. c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if(iden.eq.2.and.iii.lt.npmax) then iii=iii+1 pppx1(iii)=px1 pppy1(iii)=py1 pppz1(iii)=pz1 eee1(iii)=e1 jjj1(iii)=j1 iibb1(iii)=ib www1(iii)=dcorr*gamow wbar1=wbar1+www1(iii) pppx2(iii)=px2 pppy2(iii)=py2 pppz2(iii)=pz2 eee2(iii)=e2 jjj2(iii)=j2 iibb2(iii)=ib www2(iii)=dcorr*gamow wbar2=wbar2+www2(iii) endif c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c Edit smear to account for detector resolution. call smear(px1,py1,pz1,psx1,psy1,psz1) e1=dsqrt(mass1**2+psx1**2+psy1**2+psz1**2) call smear(px2,py2,pz2,psx2,psy2,psz2) e2=dsqrt(mass2**2+psx2**2+psy2**2+psz2**2) vs=v_source call decompose(mass1,psx1,psy1,psz1,e1, + mass2,psx2,psy2,psz2,e2,kx,ky,kz,mmom,iconvpass,vs) call kkfind(kk,mmom,akk,nkmax) if(kk.gt.nkmax) go to 444 c0(kk)=c0(kk)+dcorr cia0(kk)=cia0(kk)+1.d0 c These acceptance routines could depend any way on kx,ky,kz that c one would wish. c acceptance for parallel cut (outwards) call accpar(mmom,kx,kz,ky,ia) if(ia.eq.1) c1(kk)=c1(kk)+dcorr if(ia.eq.1) cia1(kk)=cia1(kk)+1.d0 c acceptance for in-plane cut (beam) call accin(mmom,kx,kz,ky,ia) if(ia.eq.1) c2(kk)=c2(kk)+dcorr if(ia.eq.1) cia2(kk)=cia2(kk)+1.d0 c acceptance for out-of-plane cut (sidewards) call accout(mmom,kx,kz,ky,ia) if(ia.eq.1) c3(kk)=c3(kk)+dcorr if(ia.eq.1) cia3(kk)=cia3(kk)+1.d0 c For 3-d correlation functions: call kkfind(ix3d,kx,akk,nkmax) call kkfind(iy3d,ky,akk,nkmax) call kkfind(iz3d,kz,akk,nkmax) if(ix3d.le.nkmax.and.iy3d.le.nkmax.and.iz3d.le.nkmax) then c3d(ix3d,iy3d,iz3d)=c3d(ix3d,iy3d,iz3d)+dcorr cia3d(ix3d,iy3d,iz3d)=cia3d(ix3d,iy3d,iz3d)+1.d0 endif c It should be easy to do any cut imaginable, by simply adding c a new array and making an acceptance filter. Be sure to do c the same thing in the denominator. 727 continue write(6,*) 'number of tries in numerator=',itry1 c ______________________________________________ c Now calculate the denominator c but only if ibb.ne.1 c _____________________________________________ if(ibb.eq.1.and.iden.eq.1) then call writesub(nkmax,amom,c0,cia0,c1,cia1,c2,cia2,c3,cia3) call writesub3d(nkmax,amom,c3d,cia3d) go to 999 endif c ______________________________________________ c c Some people like to choose particles for the denominator from c events with two particles in the denom., then mix particles c from different events. If this is the case, "iden" should equal c 2. If one does not require two particles to be in the event, c to use if for constructing the denominator, "iden" should c equal 1. c c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if(iden.eq.2) go to 1000 c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% write(6,*) 'Now starting on the denominator' bwtot=0.0 do 212 ib=1,nb bweight(ib)=bwtot+dfloat(np1(ib))*bw0(ib)*bcut(ib)* + (bmax(ib)**2-bmin(ib)**2) bwtot=bweight(ib) 212 continue do 213 ib=1,nb bweight(ib)=bweight(ib)/bwtot 213 continue c c ______________________________________________ itry2=0 do 827 ii=1,nmax c Pick the phase space points 844 itry2=itry2+1 c First pick the impact parameter randy=ran2() do 849 iibb=1,nb if(randy.lt.bweight(iibb)) then ib1=iibb go to 848 endif 849 continue 848 randy=ran2() do 847 iibb=1,nb if(randy.lt.bweight(iibb)) then ib2=iibb go to 846 endif 847 continue c Now pick the phase space points 846 j1=1+int(dfloat(np1(ib1))*ran2()) j2=1+int(dfloat(np2(ib2))*ran2()) if(ident.eq.1.and.j2.eq.j1.and.ib1.eq.ib2) go to 844 px1=ppx1(j1,ib1) py1=ppy1(j1,ib1) pz1=ppz1(j1,ib1) e1=ee1(j1,ib1) px2=ppx2(j2,ib2) py2=ppy2(j2,ib2) pz2=ppz2(j2,ib2) e2=ee2(j2,ib2) c Rotate just one point before testing. phi=2*pi*ran2() call ranrot(px1,py1,x1,y1,phi) c p1dotp2=e2*e1-px2*px1-py2*py1-pz2*pz1 mom2=.25d0*(mass1**2+mass2**2-2.*p1dotp2) if(msame.eq.0) then ptot2=mass1**2+mass2**2+2.d0*p1dotp2 mom2=mom2-.25d0*(pdotk**2)/ptot2 endif if(mom2.lt.mtest) go to 844 if(mom2.gt.0.d0) write(6,*) 'B: -mom**2>0, =',mom2 if(mom2.gt.0.d0) go to 844 mom=dsqrt(dabs(mom2)) phi=2.*pi*ran2() c Now randomly rotate the whole set. Note that if 'filter2' c has no phi dependence, these two rotations can be deleted. call ranrot(px1,py1,x1,y1,phi) call ranrot(px2,py2,x2,y2,phi) call filter2(mass1,px1,py1,pz1,e1, + mass2,px2,py2,pz2,e2,ifilter) if(ifilter.eq.0) go to 844 c The space-time coordinates are set to zero as they are not used c in the denominator. x1=0.d0 y1=0.d0 z1=0.d0 t1=0.d0 x2=0.d0 y2=0.d0 z2=0.d0 t2=0.d0 call transform(mass1,px1,py1,pz1,e1,x1,y1,z1,t1, + mass2,px2,py2,pz2,e2,x2,y2,z2,t2,mom,r,kdotr) c call smear(px1,py1,pz1,psx1,psy1,psz1) e1=dsqrt(mass1**2+psx1**2+psy1**2+psz1**2) call smear(px2,py2,pz2,psx2,psy2,psz2) e2=dsqrt(mass2**2+psx2**2+psy2**2+psz2**2) call decompose(mass1,psx1,psy1,psz1,e1, + mass2,psx2,psy2,psz2,e2,kx,ky,kz,mmom,iconvpass,vs) call kkfind(kk,mmom,akk,nkmax) if(kk.gt.nkmax) go to 844 d0(kk)=d0(kk)+1.d0 c acceptance for the parallel (outwards) cut call accpar(mmom,kx,kz,ky,ia) if(ia.eq.1) d1(kk)=d1(kk)+1.d0 c acceptance for the in-plane (beam) cut call accin(mmom,kx,kz,ky,ia) if(ia.eq.1) d2(kk)=d2(kk)+1.d0 c acceptance for the out-of-plane (sidewards) cut call accout(mmom,kx,kz,ky,ia) if(ia.eq.1) d3(kk)=d3(kk)+1.d0 c For 3-d correlation functions: call kkfind(ix3d,kx,akk,nkmax) call kkfind(iy3d,ky,akk,nkmax) call kkfind(iz3d,kz,akk,nkmax) if(ix3d.le.nkmax.and.iy3d.le.nkmax.and.iz3d.le.nkmax) then d3d(ix3d,iy3d,iz3d)=d3d(ix3d,iy3d,iz3d)+1.d0 endif 827 continue write(6,*) 'number of tries in denominator=',itry2 c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% go to 1001 c This next part of the code is used to calculate the c denominator if you are using only two-particle events. 1000 continue c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% do 927 ii=1,nmax 944 j1=1+int(dfloat(iii)*ran2()) j2=1+int(dfloat(iii)*ran2()) if(j2.eq.j1) go to 944 c This next line was added after Max Potekhin pointed out to me that c sometimes the same phase space point was the originator of c both the j1 and j2 phase space point. if(iibb1(j1).eq.iibb2(j2).and.jjj1(j1).eq.jjj2(j2)) go to 944 px1=pppx1(j1) py1=pppy1(j1) pz1=pppz1(j1) e1=eee1(j1) w1=www1(j1)*dfloat(iii)/wbar1 px2=pppx2(j2) py2=pppy2(j2) pz2=pppz2(j2) e2=eee2(j2) w2=www2(j2)*dfloat(iii)/wbar2 p1dotp2=e2*e1-px2*px1-py2*py1-pz2*pz1 mom2=.25d0*(mass1**2+mass2**2-2.d0*p1dotp2) if(msame.eq.0) then ptot2=mass1**2+mass2**2+2.d0*p1dotp2 mom2=mom2-.25d0*(pdotk**2)/ptot2 endif if(mom2.lt.mtest) go to 944 if(mom2.gt.0.d0) write(6,*) 'C: -mom**2>0, =',mom2 mom=dsqrt(dabs(mom2)) call transform(mass1,px1,py1,pz1,e1,x1,y1,z1,t1, + mass2,px2,py2,pz2,e2,x2,y2,z2,t2,mom,r,kdotr) c call smear(px1,py1,pz1,psx1,psy1,psz1) e1=dsqrt(mass1**2+psx1**2+psy1**2+psz1**2) call smear(px2,py2,pz2,psx2,psy2,psz2) e2=dsqrt(mass2**2+psx2**2+psy2**2+psz2**2) call decompose(mass1,psx1,psy1,psz1,e1, + mass2,psx2,psy2,psz2,e2,kx,ky,kz,mmom,iconvpass,vs) call kkfind(kk,mmom,akk,nkmax) if(kk.gt.nkmax) go to 944 d0(kk)=d0(kk)+w1*w2 c acceptance for the parallel (outwards) cut call accpar(mmom,kx,kz,ky,ia) if(ia.eq.1) d1(kk)=d1(kk)+w1*w2 c acceptance for the in-plane (beam) cut call accin(mmom,kx,kz,ky,ia) if(ia.eq.1) d2(kk)=d2(kk)+w1*w2 c acceptance for the out-of-plane (sidewards) cut call accout(mmom,kx,kz,ky,ia) if(ia.eq.1) d3(kk)=d3(kk)+w1*w2 927 continue c c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1001 continue c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c *********************************************** c Now take the ratio of the numerator to the denominator to c get the correlation function c c Note unit 66 is used for the output. call writesub(nkmax,amom,c0,d0,c1,d1,c2,d2,c3,d3) call writesub3d(nkmax,amom,c3d,d3d) 999 stop end