CRAB - VERSION 2.02 Oct. 21, 1996 Changes in version 2.02: Error in treatment of Correlation functions for non-identical particles corrected. Phase space points were not correctly read in before, when there were more than one impact parameter, and the particles were not identical. Good Luck working with the programs Please don't hesitate to call me at (517) 333-6438 or e-mail: pratt@nscl.msu.edu ________________________________________ DIRECTIONS FOR CORRELATION PROGRAMS ________________________________________ I. How to write dynamical-code results. II. Running Initialization code, CORRSTARTER III. Running main program to calc. correlation programs. ________________________________________ I. Your dynamical code should calculate the final phase space points of the relevant particles. If you are to do correlations where either the particles are identical or whether they are represented by the same set of phase space points, then you should write the points into a file, 'phasespace1.dat, in the following manner. open(21,form='formatted',status='unknown',file='phasespace1.dat') write(21,*) px,py,pz,x,y,z,t The coordinates x,y,z and t should reflect the position of the particle at its last interaction which occured at time 't'. By last interaction, I mean either its last randomizing collision or its last contact with a significant mean field. (Scattering off a mean field can also be considered randomizing.) Note that the particle never leaves contact with the long-range coulomb field of the residual system. This is the one area where physics is being swept under the rug. If the two particle types are different, e.g. you are doing K+ proton correlations, etc., then write phase space points for the first particle to 'phasespace1.dat' and the second particle to 'phasespace2.dat'. Also, write other information to 'phasedescription.dat' in an unformatted style. First write the number of impact parameters that are described by the phase space points. open(20,form='formatted',status='unknown',file='phasedescription.dat') write(20,*) nb Now for each impact parameter write the following information: do 123 ib=1,nb write(20,*) b(ib),bmin(ib),bmax(ib),npts1(ib),npts2(ib), + bw0(ib),bcut(ib) 123 continue 'b(ib)' is the impact parameter which describes the collision. 'bmin(ib)' and 'bmax(ib)' represent the range of impact parameters that are described by the 'ib'th set of phase space points. This is needed to do the weighting correctly. 'npts1(ib)' and 'npts2(ib)' are the number of phase space points that are written to 'phasespace1.dat' and 'phasespace2.dat' respectively, describing the 'ib'th set of phase space points. 'bw0(ib)' is any additional correction factor which may be needed to insure the proper weighting. Choose 'bw0(ib)' such that 'npts1(ib)*bw0(ib)' is proportional to the multiplicity for a single collision. 'bcut(ib)' is the factor which takes into account whether an event of impact parameter passes some sort of filter, (e.g. a central impact parameter trigger). If the chance that such an event of impact parameter 'b(ib)' will be recorded into the sample is 50%, then 'bcut(ib)' should be .5. If there is only one type of particle to be recorded, then type whatever one wants for 'npts2(ib)', the main program will neglect it. For example, if one runs a simulation 5000 times at 6 impact parameters of .5,1.5,2.5,3.5,4.5 and 5.5 fm, and the chances of passing the impact parameter filter are 100%, 95%, 85%,60%,30% and 5% respectively, 'phasedescription.dat' would look like this 6 .5 0.0 1.0 n1 n1 1.0 1.0 1.5 1.0 2.0 n2 n2 1.0 .95 2.5 2.0 3.0 n3 n3 1.0 .85 3.5 3.0 4.0 n4 n4 1.0 .60 4.5 4.0 5.0 n5 n5 1.0 .30 5.5 5.0 6.0 n6 n6 1.0 .05 If one performed the same calculaton but did not use the same number of events for each of the six impact parameters, but instead used 6000, 5000, 4000,3000,2000,1000 events respectively for b=.5 to b=5.5 fm, then 'phasedescription.dat' would become: 6 .5 0.0 1.0 n1 n1 1.0/6000 1.0 1.5 1.0 2.0 n2 n2 1.0/5000 .95 2.5 2.0 3.0 n3 n3 1.0/4000 .85 3.5 3.0 4.0 n4 n4 1.0/3000 .60 4.5 4.0 5.0 n5 n5 1.0/2000 .30 5.5 5.0 6.0 n6 n6 1.0/1000 .05 The most efficent way to run the program is to run a simulation and pick the number of events at each impact parameter proportional to ((bmax-bmin)**2)*bcut. This is only a statement of numerical efficiency, but any choice will be processed correctly by the program given 'phasedescription.dat' is correct. The files GAUSSMAKER.F and BJMAKER.F generate the appropriat files for gaussian sources and relativistic Bjorken geometries. For directions, see the comments in the files. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% II. Running the Initialization Program CORRSTARTER I attempted to make this program work in the easiest possible way. The program sets up the momentum mesh used in the main program and also calculates the necessary corrections to various wave-functions corresponding to the needed partial waves. This program also initializes quantities such as the masses of the two particles in the correlation programs and the fraction of the spin-states requiring symmetric or anti-symmetric wave functions. Currently the program is set up to handle strong-interaction corrections for only two partial waves at most. Since, this is set up for correlations at small relative momentum, it is unlikely a large set of partial waves will be required. A detailed list of what is initialized or calculated by CORRSTARTER follows the instructions. First compile the program. It does need any additional files linked to it. Second, run the program. The program will prompt for an i.d. specifying the pair of particles to be studied. The i.d.s currently set up are listed in the prompt. For example, if you are to do pp correlations, type 4. The program will also prompt for nkmax and delmom, which sets up the momentum values at which the correlation function will be calculated. If too large a value of 'nkmax' is chosen, information will overflow the arrays in the main program later on. This can be fixed easily in the main program(change the parameter 'nka', but may run into memory constraints at some point. One does not need to use a fixed step size delmom between values of the momentum. The main program will work with whatever momentum mesh you choose. You can go into CORRSTARTER and adjust the mesh so that it is finer when looking at a more structured part of the correlation functions, e.g. near 20 MeV/c for the p-p case. I do not believe that choosing nkmax > 50 is ever necessary. Usually, 20 is sufficient. The program writes the following information to unit 25 which is to be read by the main program. 'mass1' and 'mass2' are the masses of the two particles to be considered. 'wsym1','wsym2','wanti1' and 'wanti2' are the weights for symmetric and anti-symmetric pieces of the correlation function. These weights sum up to unity. For instance for deuteron-deuteron correlations, wsym1=5/9,wanti1=3/9, and wsym2=1/9. These correspond to the three spin-states J=2,1,0 which result in coupling the two spin one particles. 'iq' represents the product of the charges of the two particles. For p-p, K- - K- or pi+-pi+ correlations, iq=1. For p-n or n-n or pi0-pi0 correlations, iq=0. 'npart' is the number of partial waves which will be corrected for strong interactions. The program can only handle npart=0,1,2. If you wish to do more, help yourself. It takes some digging into the program to change it. 'l1' and 'l2' are the orbital angular momenta of the two partial waves. The program can only handle s and p wave corrections, so l1,l2 =0,1 only. If npart=0 then l1 and l2 will be neglected, and if npart=1, l2 will be neglected. Note that even if they are neglected, they must be written to unit 25 so that the format will be correct. 'nkmax' as mentioned above is the number of momentum values to be used. The array 'amom(kk)' describes all the different momentum values. The wave functions calculated in CORRSTARTER are the CORRECTIONS to the specific partial waves. If the strong-interaction potential =0, the corrections are zero. The program also writes the phase shift for that energy. This is needed for when the main prog. must calc. the wave function at large distances, and will use an asymptotic approximation to the partial wave. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% III. Running the main program. This program must be linked with COULSUBS.F and CSUBS.F. COULSUBS.F is a collection of routines used mainly for coulomb wave functions. These are my own routines and are not well described in the code. The method used to calculate the coul. wave functions is chosen to be numerically quite efficient. Since wave functions are calculated for specific values of the relative momentum, parameters or coefficients that depend only on the momentum are stored in arrays that make up the common block COULJUNK. CSUBS.F is a collection of routines used by the main program. Here is a descriptive list: 1. CORR(kk,mom,kdotr,r) This calculates the wave function squared as a function of 'kk' which gives the momentum value, mom which is the momentum, 'r' which is the relative coordinate, and 'kdotr' which is the relative momentum dotted into the relative position. 2. FILTER1(mass1,px,py,pz,ifilter) This returns ifilter =0 if the phase space point can not appear in the measuring apparatus, regardless of any azimuthal rotation. If the particle can be measured, a value of unity is returned. 3. FILTER2(m1,px1,py1,pz1,e1,m2,px2,py2,pz2,e2,ifilter) This returns ifilter=1 if both particles get measured and ifilter=0 if at least one particle misses the apparatus. 4. The routine TRANSFORM will calculate the values needed by corr, such as mom,r,kdotr along with the directional projections of the momentum, kpar,kin and kout. The momentum used is the invariant momentum. source. 5. SMEAR takes all momenta coordinates and smears them according to your detector resolution. The user should modify these appropriately. 6. The routine DECOMPOSE calculate the components of the relative momenta and the vector sum used for determining the correlation function. There are five options which can are determined via the parameter idirconv in 'cparams.f'. They are: Choose idirconv=1 for convention where kx,ky,kz is defined first in frame where pz1+pz2=0, then boosted outwards to the frame of the pair. In this conv. k^2=Qinv, and the z direction is parallel to the beam axis. Choose idirconv=2 for conv. where one boosts to frame of the pair, but not the outwards boost. In this conv. the z-direction is parallel to the beam axis and k^2 does not equal Qinv. Choose idirconv=3 for conv. where one boosts along the z-axis to a frame specified by vs=constant, WHICH MUST BE SET BY HAND IN THE SUBROUTINE 'decompose'. In this conv. the z-direction is along the beam axis and k^2 does not equal Qinv. Choose idirconv=4 for conv. where one boosts along the z-axis to a frame specified by vs=constant, and chooses the x-axis parallel to the total momentum In this conv. the z-axis is not necessarily parallel to the beam axis, and k^2 does not equal Qinv^2. Choose idirconv=5 for conv. where one does the same as in convention 4, only one does one more boost parallel to ptot, such that k^2=Qinv^2. This routine spits out kx,ky,kz. ky is always perpendicular to both the total momentum and the beam axis. 'mmom' is the vector sum of the three components, and it is assumed that this is what is used for binning the rel. momenta for creating correlation functions and that kx,ky,kz are used only for cutting. One can of course change this procedure. 7. If cuts are made regarding the direction of the relative momentum, these are accomplished by the routines, ACCPAR, ACCIN and ACCOUT, which make cuts for k parallel to K, k in plane and out-of plane respectively. The reaction plane is defined by the total momentum of the pair and the beam axis. 8. RANROT rotates phase space points about the beam axis. 9. KKFIND finds the appropriate value of kk, or the appropriate value of amom(kk) given the relative momentum of the pair. 10.RMARIN initializes the random # generator RANMAR and GAUSS generates gaussian numbers according to the distribution exp{-.5*x^2}. 11.CHEAPSET set the parameters in the block CHEAPJUNK which are used only when the wave function squared is calculated in Breit-Wigner parameterization which is simple. Be sure to set the parameter 'icheap' in the file 'params.f' to 1 if your are doing such a calculation. For d-alpha, K-pi or K+K- this method is used. Although the concept of a wave function is not valid in these examples, the parameters can be chosen to be consistent with assumptions of chemical equilibrium. The MAIN program is pretty well commented. Choose the parameters in 'cparams.f' 'nka' 'nba' and 'nphase' to be as small as possible without being less than 'nkmax', the number of momentum values, 'nb' the number of impact parameters studied, and 'np1(ib)' or 'np2(ib)' which are the number of phase space points which pass the single-particle filter FILTER1 which is described above. The parameter 'ident' is set to 1 if both particles are represented by the phase space points on 'phasespace1.dat', and is set to 0 if the particles are different such that you use both 'phasespace1.dat' and 'phasespace2.dat'. Be careful, about whether you want to divide out impact parameter averaging effects. If there is structure to the impact-parameter averaging then you should not divide cN(kk) by ciaN(kk). Note that doing the division reduces the statistical fluctuation in the correlation function since it compensates for the numerator and denominator having different number of points. Note on Normalization: The normalization is done such that there are the same number of accepted pairs for both the numerator and denominator. You may wish to adjust this according to your desires. The program has also been recently modified for those groups who have the desire to construct the denominator by choosing singles from events where at least two particles have been detected. To do this choose the value of "iden" to be equal to 2 instead of 1. This is set in a parameter statement near the beginning of the main program. Both the NA44 group and E859 do this. This changes the definition of the correlation function somewhat by biasing the singles towards tracks which could positively interfere with other tracks. I do not recommend constructing correlation functions this way, as it is a theoretical complication. If one wishes to Gamow correct the results set the parameter "igamow" to 1. If one does not want to do this set it equal to zero. Good luck and don't hesitate to contact me for assistance. Scott Pratt National Superconducting Cyclotron Laboratory Michigan State University East Lansing, MI 48824 PS - I would like to thank John Sullivan, Richard Morse, Mike Lisa, W.G. Gong, Nu Xu, Vince Cianciolo, Ron Soltz, Max Potekhin, and Catherine Mader for finding errors or making suggestions for improvements in the program.