TRIUMF: Canada's national laboratory for particle and nuclear physics

MONTE CARLO simulation of experiments : an introduction illustrated by the computer programs SIMUL8 (μ^+→e^+nu_e… Opat, Geoffrey I. Jul 31, 1977

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


51833-TRI-1977-04.pdf [ 16.41MB ]
JSON: 51833-1.0107830.json
JSON-LD: 51833-1.0107830-ld.json
RDF/XML (Pretty): 51833-1.0107830-rdf.xml
RDF/JSON: 51833-1.0107830-rdf.json
Turtle: 51833-1.0107830-turtle.txt
N-Triples: 51833-1.0107830-rdf-ntriples.txt
Original Record: 51833-1.0107830-source.json
Full Text

Full Text

TRIUM FMONTE CARLO SIMULATION OF EXPERIMENTSAn Introduction Illustrated by the Computer Programs SIMUL8 (y+^-e+Ve^yY in TINA and MINA Nal Detectors) and MONTE-PI ON (7r+->e+vey in TINA and MINA Nal Detectors)Geoffrey I. OpatVisitor to TRIUMF and Physics Department University of British ColumbiaOn leave from School of Physics University of MelbourneMESON F A C IL I T Y  OF:UN IVERS ITY  OF ALBERTA SIMON FRASER UN IVERS ITY  UN IVERS ITY  OF V IC TOR IAUN IVERS ITY  OF BR IT ISH  COLUMBIA T R I - 7 7 - ^TRI-77-4TRIUMFMONTE CARLO SIMULATION OF EXPERIMENTSAn Introduction Illustrated by the Computer ProgramsSIMUL8(u+->-e+VgTyY in TINA and MINA Nal Detectors) MONTE-PION(ir+-*-e+veY in TINA and MINA Nal Detectors)Geoffrey I. OpatVisitor to TRIUMF and Physics Department University of British ColumbiaOn leave from School of Physics University of MelbournePostal address:TRIUMFUniversity of British Columbia Vancouver, B.C.Canada V6T 1W5 July 1977ABSTRACTAn outline of the technique of computer simulation of particie physics experiments by the MONTE-CARLO method is presented. Useful special purpose subprograms are listed and described. At each stage the discussion is made concrete by direct reference to the programs SIMUL8 and its variant MONTE-PION, written to assist in the analysis of the radiative decay experiments u+-*e+ve\71JY and Tr+-+e+vey, respectively. These experiments were based on the use of two large sodium iodide crystals, TINA and MINA, as e and M detectors. Instructions for the use of SIMUL8 and MONTE- PION are given.INTRODUCTIONThe following reactions have been, or will be, investigated at TRIUMF using the large sodium iodide detectors TINA and MINA:-ir p -> O n-YY+ O +MM -*• it e vir+ -*■ e+vy 3y+ -*• e+vVy 4U+ -*■ e+y 5All these reactions result in two energetic particles; either two gamma rays, or an electron and a gamma ray. With the initial particle being at rest, the final configu­rations of interest have the momenta of the two most energetic particles almost anti­parallel .Reactions 2-5 are accompanied by normal muon decaysy+ -*• e+vv 6Reaction 6, having a single energetic detectable particle is recorded as a "singles" event. Reactions 1-5 having two such particles are recorded as "doubles" events.Consider for example reactions 4 and 5 which were investigated with the apparatus shown in figure 1. Scintillation counters are denoted by the symbols #1, #2, ..., #10. Collimator apertures are denoted by #20, ..., #23. A degraded (and hence multiply-scattered) pion beam is stopped by counter #3, the stop signal being 1*2*3’4. The y+ is produced by the decay ir+ y+v which gives a delayed pulse in #3. An electron is recorded as having entered TINA by the logic signal 5-7-9-(TINA). A gamma ray is recorded as having entered MINA by the signal 6*8’10*(MINA). Doubles events with the roles of TINA and MINA interchanged have a suitably modified logic.Whereas the logic signals assist in identifying the events and loosely restrict­ing their geometric configuration, the pulse-height analogue signals from TINA and MINA are the only truly kinematic variables recorded. These pulse-heights are related to the electron and gamma ray energies Ee and respectively, by the response functions ofTINA and MINA. A typical response function is sketched in figure 2. (See Appendix I for a detailed discussion.)23Fig. 2. A calibrated Sodium Iodide detector response function. R(E,E0)dE is the probability that a particle of energy EQ will give a pulse height corresponding to the range (E, E+dE).The angle between the electron and gamma ray directions 6eY is kinematically restricted by the requirement that the invariant mass of the undetected neutrino pair be positive. SpecificallyMax CO, (x+y-l)/xy) <  Sin^e^y) < 1whereEp/E,e ' Lmax 8aY = E /E '  X max 8bmax = % %The kinematic situation is shown in figure 3.8cFig. 3. The kinematically allowed region of (x,y) for a given The labels on the curves are 0gY.The allowed region for 9g = 1400 is shaded in.Some authors interchange these definitions of x and y.4Events from the putative decay u+-»e+Y (reaction 5) only populate the point (x,y) = (1,1). Events from radiative muon decay (reaction 4) are dense near the point (0,0) and along the line y = 0, and are sparse near the point (1,1). In fact the number of events in the interval l-e<x,y<l is proportional to e6. In addition, however, the finite widths of the TINA and MINA energy response functions "smear" the x-y plane distributions significantly, and thus limit the separability of events from reactions 4 and 5.If the objectives of the experiment outlined above are (i) to measure the y+-*e+Y branching ratio and (ii) to improve our knowledge of the radiative muon decay (branching ratio and n' parameter), what information would an experimenter want to have to assist in the planning and data analysis stages of the work?A suggested list would include:(a) Singles detection probability per stop, i.e. effective solidangles.(b) Singles energy distribution, i.e. a Michel spectrum convolutedwith the energy response function.(c) Doubles detection probability for each region of the (x,y) plane.(As can be seen from figure 3, each point in the x-y plane is associated with a range of angles 0eY. Given the extended target (#3) and the disposition of the counters (#7, #8) and the apertures (#20, ..., #23) the degree of acceptance of events into TINA and MINA is far from simple to compute. This is true even close to the point (1,1) where the e and y are antiparallel).(d) The x distribution and y distribution of the (cut) x-y plane.(e) The sensitivity of the above distributions to geometricmisalignments of the apparatus.(f) The sensitivity of the above distributions to finite energyresolution.(g) The sensitivity of the above distributions to beam distribution onthe target.(h) Assessment of background rates from such reactions as ir+->-e+vy, TT+-*"ir0e+v^. (A good assessment helps in the recognition of backgrounds.)This is what the programs SIMUL8 and MONTE-PION are about. However, first we shall need a little theory.II. Event Detection Probability Formulated MathematicallyWe consider here the detection of radiative muon decay events (reaction 4) by the apparatus of figure 1. We seek to compute the probability, w, that we detect the5radiative decay of a single muon, given a set of "cuts" on the energy variables x,y.Let T be the total muon decay rate. Let dTy be the differential radiative decay rate. The differential branching ratio, db, is given bydTy ^ dfle dn^db = ---  = f(e,y) • ---  • ---  • dx • dy 9N TOO 4irf(e,y) is a function of the electron and gamma ray momenta e and y given by some theoretical expression, (function subprogram THEORY in fact). db is a differential of the 6th order.Let dPg be the differential probability that a beam particle will emerge from a certain point of the degrader into a given direction and penetrate a thickness t of#3 before stopping in #3. (dPg is a differential element of the 5th order. The distri­bution associated with this element is generated by BEAMIN.)Let dPg^ be the probability that the counter i gives out pulse height energy x' when an energy x is deposited.i = TINA or MINA, or T and Mi.e. dP£i = R^(x,x')dx' 10Let a^ (e) be an aperture function for the aperturea^(e) = 1  if e goes through #j= 0  if e does not go through #jLet C(i,x',y') be a cut function corresponding to variables i = T or M;x',y' lying in some specified rangeC = 1 if variables are inside cut12= C if variables are outside cut With the aid of figure 1 we write for ww = Jdb ,dPB{a7(e)a21(e-)a20(Y-)a22(Y-)RT (x,x')RM (yjyl)C(T'xlyl)dx'dyl +13a8 Ce) a20 (e  ^a21 Cy:) a23 Cy:) SiCx *x ' ^ RT Cy»y '5 C CM*x ’y ') dx ’dy ’}We note in passing that eq. 12 contains a 13-dimensional integration with quite a complex boundary.6III.In essence (but not detail) eq. 13 has the formw = J*f(x)dp 14where dp is a multidimensional element of probability. We turn to the investigation of the evaluation of such forms by the MONTE CARLO method.General Theory of the Monte-Carlo MethodConsider the integraly = J*f (x) P (x) dnx = Jf (x) dp 15In this expression f(x) is a given function called the weight function. P(x) is a given function with the propertiesP(x) :> 0 16a/p(x)dnx = 1  16bWe can think of P(x) in several different ways:i) dp = P(x)dnx represents the probability that a point is to befound in the range dnx at x. In the present work this is how we shall interpret the function P(x).ii) P(x) is the Jacobian of a co-ordinate transformation.iii) P(x) is just simply a positive, normalised fixed function.iv) Sometimes P(x) = V y  = constant, V = volume of integration.We suppose that we have at our disposal (as a computer program) a random point generator which generates the n-tuples (xj,x2,..,xn) such that the probability dp of finding a point in the range dnx at x is P(x)dnx.We consider two statistics Sj, S2 formed from N points drawn from the random point generator. Explicitly1 N +1=1- S i 2 |177Before estimating the values of Si and S2 we define the mean and variance of the function f(x) by:y i Jf(x)P(x)dnx 19aa2 = f (f(x) -y)2P(x)dnx = Jf2 (x)P(x)dnx - y2 19bConsider the expectation valuesl N<S 1> = < f (xi)> = < f(x)> = /f(x)P(x)d2x = y 20i=li N<(Si-y)2> = <Sj2> - y2 = —  <f(xi)f (xi)> - V*2N2 i,j=l 1 J= — | N < f 2 (x)> + N(N-l) <f(x) > 2 l - y2 = 02/N 21N2 ' /Thus the statistic Si has mean y, and standard deviation a /A i . Equation 17 exhibits Si as the sum of a (large) number of independent terms, and therefore (by the central limit theorem) we conclude that Si has a normal distribution (of known mean and variance).Notes i) If N is large and a is small it is highly improbable that S 1will deviate from y appreciably. We therefore regard Si as the computed value of y with "error", ± o/Ai.ii) In general we do not know a for the function f(x) and thuswe do not know how well we have determined y. We now turn to statistic S2 for assistance with this problem.Consider the expectation value of statistic S2 of eq. 18b.E CARL S NO I | E < f2Cxi)> - <Sl2>|= ^  ju2+a2- <Si2>|= a2/N 22The last equality comes from equation 21 in which it is shown (inter alia).<Si2> - y2 = a2/NSummary Si ± v^ Sj is regarded as the computed value and error of y, in the sense of a normal distribution.Note 1 As the results show, the computed accuracy of the integral of eq. 1 does not depend on the dimensionality, n, of the space of integration. This fact underlies the reason for the importance of Monte-Carlo methods in multi-dimensional problems.Note 2 A computer program looks like:Notes: Cl) SUM1SUM2 .,, , , - . .All double precisionS2(2) Loop time must be kept to a minimumFig. 4Technical MattersWe now turn to four problems of MONTE-CARLO programs. These are, random distribution generators, specification of integration boundaries, cuts or subintegrals, and variance reduction.(i) Random Distribution GeneratorsThe most important random number generator is a computer function which generates a distribution uniform in the interval (0,1). i.e. if is the value of the function on the n ^  call, 0<x*:1.0, and xfi has probability dx of falling in the interval (x,x+dx). The successive values xn >xn+i are correlated. The FORTRAN compiler at theUBC IBM installation recognizes the statement:X = FRAND(0*0)The argument of FRAND is dummy. The sequence jXn | is determinate, but can be varied by one initial callZ = RAND(Y)The value of the real, Y, determines the sequence that will be produced by FRAND.If a different sequence jXn J is required on successive computer runs, the value of Y should be changed. The FRAND call time seems to be about 9 ysec.Although other random number generators are often supplied by the computer library, and are useful on occasion, the main programming task is to transform the uniform distribution in (0,1), into the multidimensional distribution P(x).ii) Integration BoundariesIn a multidimensional space, the region of integration may be easy to specify verbally, but impossibly difficult and/or time-consuming to specify analytically. For example (with reference to figure 1), the set of all event co-ordinates corresponding to a beam particle stopping in #3 with an electron traversing #21 and #7 and a gamma ray traversing #20 and #22, constitutes a point within the region of integration. The edge of the (13-dimensional!) region so specified defies analytic description.(Actually, another 13-dimensional region has to be added when the particles go into the opposing counters.)The way around this problem is to replace the original problem with a relatedone.Originally we had (in a simplified notation)y =ff(x)P (x)dx 23a'r°2 = J f 2 (x)Pr (x)dx - y2 23brwhere f pr (x)dx = 1 23cThe region of integration is denoted by r. It is the one with an "evil" boundary.Let R be a region with a "good" boundary such that r C  R. Within R we extend the probability function to UF XP I in a way which satisfies10andpr (x) = ^  nF XPI , x e r c  R 24bThe normalization constant, h, is the probability that a point of R lies within r.We call h the "hit" probability i.e. chance a point of R hits r.h = J P R (x)dx 25We introduce also the "miss" probability m.m = 1-h 26Further, we introduce the boundary function Br (x).B (x) = 1 x e r 27= 0 x t  rWriting programs for such functions is one of the programmers tasks, it is usually nottoo difficult, as Br (x) can usually be represented as the product of several simple factors.We introduce the new weight functionF(x) = f(x)Br (x)/h 28The new problem has mean M and variance E2,M = f F(x)PR (x)dxR 29= y (as before)E2 = f XtXP IroI2P (x)dxJR 30= (a2+my2)/h (worse than before)If h is small the increase in error may be intolerable unless N is increased signifi­cantly.Problem The programmer has to decide whether it is better to waste the machine's time by doing a large N run, or his time by slaving over an analytic boundary specification (which may also be costly to compute). Whichever way things go, you cannot just go on increasing N. Sometimes the solution comes by hybridising the problem, in some directions in R-space specifying the boundary analytically, in other directions extending the boundary a little.11Shrinking Quite often it is the reverse problem which has to be tackled. Suppose the weight ffx) contains a factor b(x) which goes to zero in some domain, then one would want to shrink the region r down onto this domain to improve efficiency.Case Studies from SIMUL8(a) Analytic specification and shrinkingXi % x '< x 1 u 31a31bMax(0,(x+y-l)/xy) = z \  * z < 1 31cwhere z = Sin2Js0ey 31dLet f(x,y,z) be zero whenever x,y,z lie outside the ranges specified byequations 31a, b, c. Equations 31 a, b, c define the region r.321 J. 1r dx J dY f dZ F(XYZ)4  0 *'o 33where F(xyz) = f(xyz)(xu-x1)*(yu-y1)*(l-z1(xy)). 34With o * X,Y, Z * 1and the relations 35 below we have analytically shrunk the spaceonto the space r.x = xj + (x^x^X 35ay = yi + (yu-yi)Y 35bz = z^x.y) + (l-zxfx.y^Z 35cWe need to generate a compatible set of variables x.y.Cose^ such thate + K y  y u1 1 1  1 1 1  y = J d x ^  dyj£  ey f(xyz) = f dx dy J* dz f(xyz)jo '<  x,y ^ 1; -1 3 Cos0eY ^ l| = R12(b) HybridisationThe integral over electron and gamma ray directions occurs as a subpart of the program SIMUL8./dfi fdO, j.,* ... , Y I _e f(S,?) 364ir J. 4tt 4TT 4ttWhere f(ft,Y) is defined as zero when (see figure 1) the electron does not traverse #7 and #21, and the gamma ray does not traverse #20 and #22. Obviously we are wasting points if we allow the electron direction 6 to be generated in all 4tt directions as the chance of hitting #7 is small. We thus transform the integral to an integration over the area of the disk of #7.dQ * . , JA e _ e*fl A # dATOO 4irr2 AWhere r is the distance from the interaction point to the electron entrypoint on #7, A is the area of the counter and ft is its normal (which points away fromthe source of radiation). Thus the solid angle probability element d!2e/4ir has been transformed into the areal probability element dA/A. i.e. dA/A is the chance an electron which does fall on A, falls inside dA.dtiy/4TT is computed regarding the electron direction to be the "Z-axis".dft dd> dCos0  Y _ ey4ir 2iri-e. =  .  ey 3gThe probability element dCos0e^/2 has been discussed in the previous section(a), the probability element d<f> /2ir was computed with the aid of RANDIR(CS,SN) whicheygenerates the sine and cosine of a random direction in the plane.With these changes the new weight function is:F(S,?) = * £ C ^ )  39In the region of energies x,y=l the electron and gamma ray move almost antiparallel. With the electron constrained to traverse #7, we find the overall hit probability for the entire event to be h-0.8, only a small source of variance increase.13iii) Cuts and SlicesIn most applications one wants to know, not only the overall detection rate, but also the detection rate corresponding to say a particle energy in some subrange. Such events are "cut" out as they occur in the MONTE CARLO loop of figure 4 and their subtotals are also formed. In fact, a very large number of one and two dimensional cuts are often performed simultaneously within the main loop in order to display the functional dependence on the selected variables.From the point of view of a single cut it is equivalent to the replacementin eq. 15f(x) -* f(x)Cfx) 40where C(x) = 1 x e c c r= 0 d e c c r  41and c is the region of the cut within r.Analogous considerations to those given previously show:Uc = hy^ 42aa2c = h(o^2 + my^2) 42bIn eqs. 42, h and m are the cut hit and miss probabilities and y' and a ' 2 are thec cmean and variance of f within the cut. Quite often o| * 0 and m = 1. In this case the fractional uncertainty in yc is given bygc „ 1y Vii y jh ti 43which is simply the expected Poisson partitioning error. In fact eq. 43 gives a minimum value to the error within a cut.iv) Variance ReductionIt is clear that if the variance of f(x) is large the MONTE CARLO method will work poorly. Several devices are available to reduce the variance and this subject forms much of the content of the literature on the method (Refs. 1, 3). One device is simple if it can be applied.14Suppose a function Q(x) can be found, such that f(x)/Q(x) has a small variance. Then the replacementf(x) -*• f(x)/Q(x) = f ' (x) 44aP(x) -> Q(x)P(x) = P' (x) 44bcan often give a great improvement. The present program chains do not have explicit variance reduction included. This is partly because of the requirement of uniform (fractional) accuracy within each cut as well as good overall accuracy, and partly because of its difficulty of implementation.One dimensional example of problems (i) - (iv)We compute the trivial integral 1f 2xdx = 1 J0using the MONTE CARLO method.100 Ei=l, UUS 1 =  1 0 0  E 2 x i With 0 '< x ± '< 1This process was repeated 1000 times using a different set of random numbers each time. The distribution of values of Si is shown in figure 5a. The mean and variance of the 1000 numbers so found is in good agreement with that expected from the theoretical evaluation of the variance.In figure 5b the effect of extension of the boundary are manifest. Note to compute this case the following sum is found,. 100s l = I5o E 2 ( 2xi)0(1-Xi). with 0 * Xi '< 2 i=lIn figure 5c the variance reduction technique is illustrated. The variable was y = Coshx e (1, Cosh(l)).In figure 5d the larger variance within a cut is illustrated.Figure 5.16IV. The Programs SIMUL8 and MONTE-PION(i) General Comments (see figure 6)The programs are designed to be efficient within the main loop. With minimum options in effect the loop time is about 8 millisec on the UBC IBM system, and rises to 15 millisec with all options in effect. A typical run generates 10,000 (NSTAT) statistics. This number of statistics produces about 1% or better accuracy in the "globally" averaged quantities and about 10% accuracy in small cut quantities. The limits to the electron energy are set by XMIN and XMAX and the gamma ray energy by YMIN and YMAX.* The number of subdivisions within these regions is controlled by NX and NY. If high precision results are required for any subregion it is much better to run the program with say NSTAT = 10,000 statistics all falling within the region of interest, rather than with a massive run of say NSTAT = 1,000,000 falling in a wider region.The following techniques have been used to increase the speed of the program:a) Avoiding lengthy algebraic and transcendental operations in the mainloop (where possible).b) Reducing the number of calls to subprograms (where possible).c) Writing initialisation sections within each subprogram. Whereappropriate, all constants are computed on the first entry to the subprogram. A switch, INIT, is then set to prevent subsequent calls entering this region again.It was found that by using the FORTRAN H compiler only a 10% improvement in speed was achieved. Due to the uncertainties surrounding this compiler its use is not recommended.(ii) The Co-ordinate systemThe programs assume a right handed cartesian co-ordinate system. With the exception of the subroutine BEAMIN any origin and orientation of the co-ordinate system is satisfactory.Lengths are all specified in centimeters.Angles are all specified in radian.Energies, masses and momenta are in MeV, MeV/c2, MeV/c respectively, or aredimensionless fractions of their maximum values.In BEAMIN (as it is presently written) the co-ordinate assumptions are:-i) Z-axis is parallel to the beam direction.ii) X-axis is the axis of TINA and points into TINA.iii) Y-axis is vertically up in the laboratory.*See the Appendix I for an elaboration of relationship between energy scales and pulse heights.Main loop17iv) The outside iron face of TINA is the Z=0 plane (see figure 1).Recommendation: In choosing a co-ordinate system do not necessarily choose themost symmetric or beautiful one, choose one that is tangible and has easy access on the experimental floor!MAIN PROGRAM SUBROUTINES USED(other than those provided by system default options)COMMON BANKSRUN (Main)INDATARUNCNTRSBEAMBEAMIN N0RML2BEAMCNTRSEVENTENADJCNTAPT( ;EVTBNK)CNTPNT( ;CNTRS,EVTBNK) RANDIRTHEORY or DGDXDYDEDXRANGERESOL(RANVAR;TINMIN)RUNCNTRSEVTBNKTINMINEVTBNKDATSUMRUNEVTBNKSUMBNKPNTOUT GRAPH RUNSUMBNKSTOPFig. 4. Structure of the Programs SIMUL8 and MONTE-PION.18iii) Running the programThe source programs are on files SIMUL8*S and MONTE-PION'S, the object programs being on SIMUL8*0 and M0NTE-PI0N*0. A typical run from BATCH has a card deck which looks like:Card 1 Service Request Card $SIG OPAT T=20 P=20 Password $RUN SIMUL8-0R... (run card data) compulsory(beam card data) _ compulsorycounter or aperture data6+NCNT7+NCNT8+NCNTB..C.. C.. C. , Any number of cards, (NCNT) in any order, bearing counter or aperture numbers in the range 1 through 30.C... $END $SIGiv) The Data Cards - their content and formatThe data cards are all read by subroutine INDATA which echoes their contents on the line printer in the same format.The RUN card has the following data string:Name NSTAT XMIN XMAX YMIN YMAX NX NY IESWeg. R 10000 0.500 1.000 0.300 0.900 10 12 1FORMAT (1X,I9,4F10.3,3I10)NSTAT = number of statistics generated*XMIN, XMAX, YMIN, YMAX set lower and upper limits to the variables X,Y.X,Y refer to electron and gamma ray energies or pulse heights.*NX, NY are the number of bins into which the X and Y distributions which appear in thegraphs and tables will be cut. 1 < NX,NY '< 20.IESW is a mode switch.For IESW = 0  No allowance is made for energy loss or the energy response of TINA or MINA. The energy adjusting subroutine ENADJ is skipped.= 1 Energy loss and sodium iodide response is allowed for.= 2 The initial undegraded values of the electron and gamma ray energiesare frozen to be equal to XMAX and YMAX and their directions anti­parallel. This option enables the monochromatic energy response to be studied. (Available in SIMUL8 only.)*See the Appendix 1 for a fuller description of the relationship between pulse heights and energies.19The BEAM card has the following data string:Name IBMTYP BM(2) BM(3) BM(4) BM(5) BM(6) BM(7) BM(8)eg. B 4 -7.620 1.270 5.000 4.445 0.077 10.000 5.000FORMAT (IX,19,7F10.3)The meaning of the items of the data string BM(2)-BM(8) depends on the type of beam (IBMTYP). See the listing of BEAMIN for a description of the beam types available and the meaning of the variables.The COUNTER cards (or aperture cards) each have the following data string:Name IC CNT(1,IC) CNT(2,IC) CNT(3,IC) CNT(4,IC) CNT(5,IC) CNT(6,IC) CNT(7,IC) CNT(8,IC)eg. C 3 -7.420 0.000 0.000 0.881 0.000 0.474 20.320 0.635FORMAT (IX,I2,F7.3,7F10.3)The index number IC is the same as that used on the experimental floor with thefollowing restriction, 1 '< IC '< 30.For the counter (or aperture) with index number IC the data string has the following meaning:CNT(1,IC), CNT(2,IC), CNT(3,IC) are the X,Y,Z co-ordinates of the point at the centre of the counter or aperture (on the entry side in the case of a thick counter). CNT(4,IC), CNT(5,IC), CNT(6,IC) are the X,Y,Z components of the unit vector normal tothe plane of the counter. The sense of this unit vector is from the entry to the exitside of the counter. (The sense of this vector is used to test whether a ray pene­trated the counter from the expected side.)CNT(7,IC) is the RADIUS (not diameter) of circular counter or is the side length of a square counter.CNT(8,IC) is the thickness of the counter or the thickness of material (if any) blocking an apertureNote: This 8 parameter specification of counters and aperture could in principlebe incomplete if:i) they were all made of different materialii) they were not circular or square in shapeiii) their orientation (for non-circular) counters required specification,v) The Named COMMON BanksThere are no unnamed COMMON banks. The contents of the named COMMON will nowbe described.C0MM0N/RUN/NSTAT,XMIN,XMAX,YMIN,YMAX,XRANGE,YRANGE,NX,NY,IESW These variables are all on the RUN card except for XRANGE = XMAX-XMIN YRANGE = YMAX-YMINCOMMON/BEAM/BM(8)EQUIVALENCE (IBMTYP,BM(1))These variables are all on the BEAM card.20COMMON/CNTRS/NCNT,CNT(8,30)These variables are all on the COUNTER or aperture cards except for NCNT, which is the total number of counters and apertures in the system. C0MM0N/EVTBNK/P(3),TB,E(3),G(3),X,Y,COSEG,WT(10),KET,KEM,KGT,KGM,THICK EVTBNK holds the variables for the current event.P Cartesian coordinates of interaction point.TB Thickness of material penetrated by beam in stopping counter beforecoming to rest.E Unit vector in electron direction.G Unit vector in gamma ray direction.X,Y Electron and gamma ray energies in units of the maximum energy,before and after energy adjustment due to slowing and resolution by ENADJ.COSEG Cosine of the angle between the electron and gamma ray.WT( ) A string of weights calculated in EVENT and used by DATSUM.WT(1) Weight for u+-»e+v_v (Michel Spectrum).6 MWT(2) A phase-space-like weightWT(3) Weight for y+-*e+v v y in SIMUL8 ) calculated using+ + . MriMTC DTnM } THEORY or DGDXDYir -*e vgy in MONTE-PION jWT(4) = GEOMWT Geometric weight to do with electron back counter solidangle; see CNTPNT.WT(5) =0,1 TINA/MINA front aperture function for electrons.WT(I) 6 * I '< 10 are free.KET,KEM =0,1 Depending on whether CNTPNT put the electron through the TINA or MINA back counter.KGT,KGM =0,1 Depending on whether CNTAPT put the gamma ray through the TINA or MINA front and back apertures.THICK Is the total thickness of material penetrated by the electronallowing for its (slant) trajectory out of the stopping counter and through the various other counters.C0MM0N/SUMBNK/A(15),AX(20,15),AY(20,15),AXY(20,20,15),NC0UNT DOUBLE PRECISION A,AX,AY,AXYA,AX,AY,AXY are the arrays in which DATSUM assembles the partial sums corresponding to the overall rate, X-distribution, Y-distribution and XY-distribution.NCOUNT is incremented by unity on each entry to DATSUM. When NC0UNT = NSTAT the arrays are processed to be ready for output by PNTOUT.21V. The Subroutines and FunctionsEach subroutine and function is described below, and where it is appropriate any special algebra, geometry or listings are displayed. The subprograms are discussed in alphabetical order for easy access. For full listings consult R. Poutissou of TRIUMF and/or M. Hasinoff of the Physics Department, U.B.C..BEAMINA call to BEAMIN returns the cartesian coordinates P(l), P(2), P(3), of abeam particle stopped in counter #3, (see figure 1), generated in accordance with anassumed beam profile probability distribution. It also finds the thickness, TB, of matter penetrated by the beam in Counter #3. The centre line of the beam is assumed to be parallel to the positive Z-axis. The type of beam profile is determined by avariable IBMTYP which is listed on the Beam card.tT f iIFig. 7. Beam stopping in counter #3.The following types of beams exist at the time of writing:IBMTYP = 1 Parallel beam, uniformly distributedwithin a rectangular cross-section.= 2 Parallel beam, uniformly distributedwithin a vertically cut-off circularaperture.= 3 Parallel beam, cut-off by a circularaperture with gaussian profiles in the horizontal and vertical directions of differing width.22IBMTYP = 4 Beam multiply scattered by a degrader into a short collimator. Beam profile on the degrader has gaussian profiles in the horizontal and vertical directions of differing width.the most complex we describe it further in the form of several notes.Beam card variables IBMTYP = BM(1) described above = 4BM(2),BM(3) = (XCENT,YCENT) coordinates of the centre line of the beam.BM(4) not used.BM(5) radius of collimator.BM(6) multiple scattering angle 80 as used in exp-(0/20e)2. BM(7),BM(8) full width at half height of the beam profiles in the X and Y directions respectively.N0RML2 is used to generate both beam position on the degrader (XX,YY) and the multiply scattered beam direction vector (XXN,YYN,1.0).The beam stopping position is chosen uniformly through the counter #3.The following equations are employed to find the intercept of a ray on a plane.r = a + X(1 (Equation to ray through a with unit direction ft). (r-rQ)*ft = 0 (Equation to plane through rQ with unit normal ft). Thus the distance, X, from a to the plane is given by X = (rQ-a) *ft/(ft*ft) and r = a + (rQ-a) «ftO/(fl*fl).CNTAPT(IC,P,U,KOUNT,DX)This subroutine tests whether a ray starting at point P moving in unit direction U penetrates counter or aperture, IC, from the correct side. If it does, K0UNT=1 and DX=thickness (in cm) of counter penetrated. If it does not, KOUNT=0, and DX=0.0.CNTPNTThis subroutine "puts an electron through" the TINA or MINA defining counter. More specifically CNTPNT:a) Decides whether to put the electron through TINA or MINA.b) Sets the flags KET and KEM accordingly.c) Finds the point on the defining counters (7 or 8) which ispenetrated.As beam 4 is a)b)c)d)23d) Computes the unit electron direction E(l),E(2),E(3).e) Computes the thickness (cm) of counter penetrated, THICK.f) Computes the geometric weight GEOMWT = WT(4), The geometric weight is the geometric factor in eq. 39 above, namely d*M/(4Trr2). The area A is the sum of the TINA and MINA defining counter areas, i.e. A = ifR2TINA + 1tr2mina- The program regards these defining counters as disjoint patches of a single common area, A.g) This subroutine (together with several others) needs two orthogonal vectors in a plane normal to a given (unit) vector. The following device is used:Let (a,b,c) be the given unit vector, a2 + b2 + c2 = 1.Defining d 5 (a2 + b2)'2 we find that (b/d,-a/d,o) and (ac/d,bc/d,-d) are two such unit vectors.DATSUMThis subroutine at each call summarises the data in a way to be described.On the last call (for which NCOUNT=NSTAT) the data arrays are normalised. DATSUM carries out the following at each call:a) Combines the computed weights, WT, with the logical variables KET, KGT, KEM, KGM, to produce the weights W relevant to the experiment at hand.b) It allocates a bin, IX, IY, to the energy variable pair x 1 and y 1. There are NX(NY) bins between XMIN and XMAX (YMIN and YMAX).c) The double precisioned arrays A(15), AX(20,15), AY(20,15),AXY(20,20,15) have added into them up to 15 different weights, not cut at all, cut into X bins, cut into Y bins, and cut into XY bins respectively.d) By dividing the theoretical weight into the geometrically modified weights, DATSUM also produces the detection probabilities for each region of x ’.y’ space.e) The letter A roughly corresponds to the statistic SI.Note: As written the statistic S2 and its cut versions are notcomputed by DATSUM, double precision arrays B(15), BX(20,15),BY(20,15), BXY(20,20,15) need adding to the SUMBNK for this purpose.DEDX(T)DEDX(T) gives the loss of kinetic energy of the electron per unit path length in a polystyrene counter. The kinetic energy T is in MeV and DEDX is in MeV/cm.The function employed is a piecewise quadratic fit by A.J. (Hannes) Barnard to the electron collision loss tables of M.J. Berger and S.M. Seltzer, National Aeronauticsand Space Administration (Washington, D.C.) report NASA SP-3012 (1964). The accuracyof the fit is adequate to the task at hand.DGDXDY(X,Y)This function evaluates the differential branching ratiodb/dxdy = (dr /dxdy)/r for the decay ir+->e+v y , and as such it is the main weight Y ®function of MONTE-PION. The formula is in accord with:i) G.I. Opat - Personal Notes 2/April/1977. (J-M. Poutissou and P. Depommier have copies).ii) P. Scheck and A. Wullschleger, Nuclear Physics B67, 504-517 (1973).iii) V. Soergel, Proc. SIN Spring School on Weak Interactions andNuclear Structure, Zuoz. Switzerland, April (1972) Vol. II, P.51.The inner bremsstrahlung term and the structure dependent terms are evaluated but the interference term is not. This omission has two origins:a) It is very small under all conditions.b) The sign of this term is not known (to me at least).ENADJThis subroutine adjusts the energy of the positron and gamma ray. For the positron ENADJ:i) Tests whether the positron has sufficient energy to penetrate all relevant counters. If not WT(5) = 0.0. Uses RANGE(E).ii) Decreases the electron energy by the collision loss term due tothe matter penetrated. Uses DEDX(E).iii) Adds the annihilation energy mgc2 to the energy.iv) Randomly re-allocates the energy in accordance with the Naldetector response function.Note: The electron energy is taken as the total energy, not as the kinetic energy.For the gamma-ray ENADJ:i) Randomly re-allocates the energy in accordance with the Naldetector response function.The resulting energies are supposed to be analogues of the pulse height put out by TINA and MINA, on a common scale, with X or Y = 1.0 corresponding to an energy EMAX.EVENTThis is a very major subroutine which works out nearly all the kinematic variables and evaluates the weights for each event. It begins with the stopped particle produced by BEAMIN and puts out its data into the bank EVTBNK.2425In order,the following sequence is followed:i) Electron energy is generated randomly within limits set.ii) Michel distribution weight is evaluated and put in WT(1).Note: The philosophy of calculating a "singles" event in parallelwith the "doubles" events is very sound, as these are also collected by the experiment. The "singles" data forms a basic check on the experiment.iii) Gamma ray energy is evaluated randomly within limits set.iv) Cosine of the angle between the electron and gamma ray is assigned.a) In SIMUL8 (y -*-e vvy) this is done randomly in accordance with the limits set by eq. 31c.b) In MONTE-PION XOO -+e vy) this angle is directly computed from the energies.v) In SIMUL8, THEORY evaluates the differential branching rate fory +e vvy. In MONTE-PION, DGDXDY evaluates the differential branch­ing rate for Tr+-» The electron (which is common to y+-*e+vv and y+-»-e+vvy or to+ + —  + + y ->e vv and tt -+e vy) is put through the defining counter(#7 and #8) by CNTPNT.vii) Aperture acceptance tests are performed on the electron.viii) The gamma ray direction is generated.ix) Aperture acceptance tests are performed on the gamma ray.Note 1 WT(1) is theoretical singles weight.WT(2) is the theoretical phase space (only) doubles weight.WT(3) is the theoretical doubles weight.WT(4) assigned by CNTPNT is the geometric weight.WT(5) is the electron front aperture acceptance weight. It is changed to zero by ENADJ if the electron fails to penetrate the required amount of matter.Note 2 In SIMUL8 (only) if the RUN card switch, IESW=2, the electron energy is setto XMAX, and the gamma energy to YMAX. The angle between the electron and gamma ray is set to 180°. Also WT(1) = W(3) = 1.0. This feature was introduced to assist with experimental diagnostics.GRAPH(A,AM,N)Plots the graph of N values of a non-negative function stored in the array A.The maximum value of the string in A is AM, which must satisfy 1.O^AM^IO.0. GRAPH also selects a pleasing height for the graph.26INDATAReads in data cards in formats already discussed.MAINSee figure 6. If IESW = 0, ENADJ is skipped.N0RML2(X,Y)Returns two independent variables X,Y each of which is normally distributed about zero with unit variance. The algorithm on which it is based is exact. In spite of the use of transcendental operations, it was found that N0RML2 was faster (by a•kfactor of 2) than two calls to the usual generator of unit normal variates based on the central limit theorem.SUBROUTINE NORML2<X,Y)C NCRML2 RETURNS TWO INDEPENDENT NORMAL (GAUSSIAN) VAR IABLES.EACHC OF UNIT VARIANCE CENTRED ON ZERO.N=0 1 N = N + 1IF(N.GT.1CCCC) GO TO 100 A=1.0—2.0*FRAND(0.0)8 * 1 . C-2.C*FRAND(0.0)RR=A**2+B**2IF((RR.GT.1.01.OR.(RR.EO.O.O)) GO TO 1 C=SQRT{-2.0+ALGG(FRAND(0.0))/RR)X=A*CY=B*CRETURN1000 FORMAT!' FAULT IN SUBROUTINE N0RML2.M0RE THAN 10000 TRIALS') CALL EXIT ENDSuch a generator function is: FUNCTION GAUSS (0.0)GAUSS = -6.0 DO 1 1=1,12 1 GAUSS = GAUSS+FRAND (0.0) RETURN END27PNTOUTOrganizes prints and graphs the contents of the SUMBNK. As presently written it is adequate for the programs SIMUL8 and MONTE-PION. By way of criticism it needs the following modifications:i) Print and graph control from a run card.ii) Better graphics, including possibly a contour plot of AXY.iii) Provision for the "error" statistic output, /S2, if this is everincluded in DATSUM.iv) Generally greater flexibility.RANDIR(CS,SN)Finds the sine and cosine of a direction in the plane, without finding the angle itself. RANDIR uses no transcendental functions at all and hence is very fast.Fig. 8. The angle a is uniform in (0,2-n).SUBP OUT INF RANOIR(CS » S N)C**»* RANDIR FINOS THE SINE AND COSINE OF A PANDCM ANGLE BY THF METHOD C CF VON NEUMANN.I A=FRAND(C.O)e=1.0-2.0*FRANO(0.0>A2=A**2e2=B**2C2=A2+B2IFIC2.GT.1.0) GO TO 1 CS=(A2-B2)/C2 SN=2.0*A*B/C 2 PETURN FND28RANGE(T)Gives the range (cm) of an electron of kinetic energy T(MeV) in polystyrene. The function RANGE was constructed in a similar manner to DEDX and from the same reference; see DEDX.RANVAR(XMIN,XMAX,N,P)This function implements a general procedure for generating a random variable for an arbitrary one dimensional probability distribution.Fig. 9. A probability distribution function P(x) is shown in (a) and its integral Q(x) in (b).Let dp be the probability of finding x in (x,x+dx). Thendp = P(x)dxdefines a probability distribution function. Let Q(x) be defined byxrQ(x) = J P(x)dx.XMIN4546It is easy to show that if y is uniformly distributed in (0,1) then x given byx = Q~1(y) 47has the required distribution.RANVAR expects P(x) to be specified at N equally spaced values of X in the interval (XMIN,XMAX). The consecutive values of P(x) are given in the array P, the 1st being P(XMIN) and the last P(XMAX).29On an initial call RANVAR integrates P(x) by passing a cubic polynomial through each set of 4 consecutive points. The end intervals are treated separately. The values of the integrals at the N points are stored in the array PI, and the differences of the integrals at successive points in the array DPI.On every call RANVAR "looks up" the value of the required random variable by linear interpolation between the points.Note: i) RANVAR could be made a little faster by using a "binary search" forthe required random variableii) RANVAR, as written, can only be initialised once, thus subsequentcalls always are evaluated using the first probability distribution function.FUNCTION RANVA°(XMIN,XMAX,N,P)C**** RANV4® RETU°NS A RANDOM VARIABLE IN THE RANGE <XMIN,XMAX) BASF D ON C A PROBABILITY DISTRIBUTION G ! V PN AT N EQUALLY SPACED POINTS.THEC RIRST PCINT CORRESPONDS TO XMIN,THE LAST TO XMAX.THE PROBABILITYC DISTRIBUTION,vJHICH NEED NOT BF NORMAL ISFD , I S TAKFN AS A C UBICALLYC SMOOTHED FUNCTION THROUGH THF POINTS.DIMENSION P(\),PI(1000),DPI<10C0)CATA INIT/O/IF( INIT.NF.O) GO TO 10 INIT = 1D P I ( 2 ) = 9 . 0 * P ( 1 ) + I 9 . 0 * P (2)-5.0*P(3)+P(A)N1=N-1C" 1 J = 3 ,N 11 DPI C J)=-P(J-2J *1 3.0*(P (J— 1) P ( J ) l-P ( J ♦ 1J DPI{N)=P(N-3)-5.0*P<N-2»«-19.0*P{N-l)*9.0*P(N)P I C 11 = 0 . CCO 2 J = 2 , N2 PIIJ)=PI(J-l)+0PI(J)CO 3 J = 2,NC P I (J 1 = C P I (J >/PIIN)3 PI(J)=PI(J)/PI(N)C X = ( X M A X - X M I N I / F L O A T ( N - l )10 R = E R A N D (C . O ) no 11 J = 2,N T = R - P I (J )I F ( T. Lc . 0. 0 ) GO TO 1211 CONTINUE12 R A N V A R = X M I N * C X * { F L O A T {J-l ) +T/DP I (J ) )RETURNEND30RESOL(E,IEORG,ITORM)The function RESOL allocates to the deposited energy E (MeV) a random variable, the pulse height RESOL (MeV), according to the pulse height response (probability) function. (See Appendix I for details of the pulse height scale.)RESOL is based on 3 hypotheses concerning the response function. If E is the deposited energy and H is the corresponding pulse height these hypotheses are:-(a) H0 E, where HQ is the most probable pulse height. (In RESOL,H0 = E).(b) AH « E*5, where AH is the full width at half height of the responsefunction.(c) The response function profile is universal.RESOL implements these hypotheses by the use of a data set stored by a BLOCK DATA statement. The data set includes:i) The four fractional widths, AH/HQ = F, for the cases electrons andgamma rays into TINA and MINA.ii) The energy EZERO at which these four fractional widths pertain.iii) The universal response function in standard form. The N(unnormalised) values of this function are held in RESFUN. The standard form requires that the most probable point have abscissa zero, and that the full width of half height be unity. The lower and upper values at which the response function goes to zero are RL and RU, (RL <0.0). See figure 10.RANVAR implements the generation of a random variable in accordance with the universal response function.'RESOLA r t  «  l ose.Fig. 10. The data representing the universal response function in standard form.31P U N C f l u N  d u C c i l u s l u c a e s l l c d bylhIC * * * *  x E b u L  S H t A X S  M t  t N t X b Y  U F  t L E L f H U N S  A N Q G A m m A  H A Y S  I n  T I N A  A N D  C "IliMA A C C U X o I  N o  Tu  I H E  t X P £ k  1 H £  n  I A l. P H U F I L E .m c o r l p S  / r I g hb 8 S ( + i -  >  U / N » X £ i i P U . N (  1 0 0  ) » F (<? * t! J V Y ( u d c - u d D P ) * u o l S  m i t G H l j = l , d  F U X  t L E L I H U N f  e D o g D r + D vC F O X  I 1 N A , M I N AC t ANij a u C G p  D P u l S  "*Ev7 = E / t 4 t H UFiF g ure I(i)* F (I tU XG ,IIORM)*HAnvah(HL»HU#Nf HESFuN)H t S U L = t ^ E X u * 2  K t I U X Nt NL)b L u C *  D A T AC U M M u N  / r i i y « l i ' j / K L » H u » N # H t b F U i y ( 1 0 u j , F l 2 # < J j # t ^ E W U # E M A X * i J E ^ l N » E A N N l i 1  U A  I A W L * H U » N , E Z t x 0 # E H A X / - a . 7 3 l 9 , 0 . 9 a b X , l 9 , 5 d . 8 d 9 / , 5 d . 8 2 9 7 /m n d u  C u S  f I) A f A i s  FH(jrt p s f u a . > p t u >  t B c >  P l O N  C E x  U A I A  C t ? b / 3 / l 9 / 7  J .O A  ) A h E S F u N / 0 . 0 , t ) . b , e » l . l , S S . o , ‘* 9 . 3 , 6 / . l , 9 < J . l , l l / . ^ , l a ^ . <»,1 d 3 3 . 9 , a y © . 9, 3 5 9 .  7, / 7 a .  9 ,  1 1 8 7 . 3 ,  1 9 1 3 . 8 ,  d 4 7 d . d , l b o a . u , 3 7 a .  5,* y  , u > b i * y . u /  M O N O  T O E E N C A R L 3i l o o b a /CC I H t  bS v A L U t S  J F  H E  S F  U N  S P t C i F Y  : p d u > 8 m D i i  v N > Y c S l , u d C 6 i C c * l c o  l c c l p uC H t S P u N S E  F U N C I l U N  U n  A S C A L E  I n  w h I C h  T H E  P E A K  I S  A I  y . O  O E S T H E  F U L L  WC w l O T H  AI H A L F  H E I G H )  I S  l . UC K L *  H U  A x E  N > 1  L U w t H  D S c  U P P E H  V A L U E S  AI w H I C H  T H E  K E S P U N S E  G O E S  T O  0 0C U N  I H I S  S L A L t  .t NOTHEORY(X,Y,DL)Evaluates the differential branching ratio for the decay u+-*e+vvy and assuch is the main weight function of SIMUL8. The formula is taken from G. KSllIn,Springer Tracts in Modem Physics 46, 96 (1968). The electron mass has been set to DL = l-Cos9eyVI. References^ M. Hammersley and D.C. Handscomb. Monte Carlo Methods (Halsted Press (Wiley) U.S.A., or Methuen U.M., 1965)NNv ) 6dw11n” * ^ -retrospective and prospective survey of the Monte Carlo Method.S.I.A.M. Rev. 12, 1-63, (1970). (This review is particularly valuable for its long list of references, some of which are to physics applications.)KleiJnen- Statistical Techniques in Simulation, Vols. I § II (Marcel Dekker, N.Y., 1974-5).F^E. James. The CERN program FOWL and its long writeup. (CERN Program Library -?heecov*r (FOWL evaluates integrals of arbitrary matrix elements overthe covariant phase space of any number of particles using MONTE-CARLO methods.)32APPENDIX 1Important Note on Energies and Pulse HeightsBefore entering ENADJ, X and Y are the total energies of the electron andgamma ray in units of the maximum theoretically allowed energy, EMAX.(EMAX = Mp/2 in SIMUL8), Thus 0.0*X,Y*1.0.In ENADJ, X and Y are converted to an MeV energy scale. The (positive)electron energy, ENE, is adjusted to allow for collision loss in the counters and the annihilation energy, whereas the gamma ray energy ENG is not adjusted at all. At this stage ENE and ENG are the energies dumped into the Sodium Iodide crystals by the electron and gamma ray.ENADJ now sends ENE and ENG through the crystal response routine RESOL. The output variables of RESOL, are also called ENE and ENG but now represent pulse heights in the standard convention.The standard convention for pulse heights is described as follows:i) Zero pulse height corresponds to zero energy dunped into the crystal.ii) The most probable pulse height when energy E is deposited is assumed to be proportional to E on physical grounds. For this reason, in the standard convention for pulse heights the most probable pulse height for energy E MeV deposited is taken to be numerically equal to E. We note that the pulse heights correspond­ing to deposited energy E can be less than or (slightly) greater than E.Finally, ENADJ divides the pulse heights ENE and ENG by EMAX to produce the dimensionless variables, again called X and Y. We note that X and Y are positive variables whose range can extend a little over unity.When IESW=0, ENADJ is not involved and X,Y remain energy variables.The run card variables XMIN, XMAX, YMIN, YMAX set limits on the ranges of the output variables X and Y.If IESW = 0 these parameters set the limits to.electron and gamma rayenergies.If IESW ^ 0 these parameters set the limits to the electron and gamma ray pulse heights. In this case the program has to compute the corresponding limits to33to the ranges of the energies. This is done in an initialisation section of EVENT and the limits are denoted XL, XU, YL, YU. These limits are chosen to be just wide enough, so that no energy point excluded could have produced a pulse height in the desired range.lam i


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items