MONTE CARLO SIMULATION OF RARE DECAYS OF T H E 7T° Reena Meijer Drees B . Sc. (Hons, co-op) University of Waterloo A THESIS T H E SUBMITTED IN PARTIAL REQUIREMENTS M A S T E R FULFILLMENT OF F O R T H E DEGREE OF OF SCIENCE in T H E FACULTY OF GRADUATE DEPARTMENT OF STUDIES PHYSICS We accept this thesis as conforming to the required standard T H E UNIVERSITY O F BRITISH COLUMBIA September 1988 (c) Reena Meijer Drees In presenting degree freely at this the available copying of department publication of in partial fulfilment University of British Columbia, for this or thesis reference thesis by this for his thesis and scholarly or for her Department DE-6 (2/88) Columbia I further purposes gain the requirements I agree shall that agree may representatives. financial permission. The University of British Vancouver, Canada study. of be It not is that the Library permission granted by understood be for allowed an advanced shall for the that without make it extensive head of my copying or my written Abstract A basic introduction to the techniques of writing computer simulations for particle physics experiments concerning rare decays is given. A very brief outline of the theoretical basis for cross-section calculations is followed by more in-depth discussions of Monte Carlo sampling techniques, including the use of the C E R N package D I V O N . The reader is introduced to the concepts involved in detector simulation; specifically, the C E R N package G E A N T in presented and discussed in some detail. Finally, methods of modelling external conversion processes are taken up. Illustrations for the discussions are taken from the various rare decay modes of the 7r°. Example code and some sample results comparing simulated and real data are given in the appendices. u Table of Contents Abstract ii Acknowledgement vi 1 Introduction 1 2 D a t a Simulation: Event Generation 3 2.1 The Differential Branching Ratio 4 2.2 Monte Carlo Sampling Methods 9 2.3 Laboratory Frame Kinematics 15, 2.4 The Event Vertex 18 2.5 Modelling External Conversions 19 3 Detector Simulation: Tracking 24 3.1 G E A N T : Introduction and General Philosophy 26 3.1.1 Detector Definition 30 3.1.2 Tracking 34 3.1.3 Digitization and Output 36 3.2 4 External Conversions and Compton Scattering Conclusions 38 43 Bibliography 45 A 47 D I V O N coded example in B G E A N T detector definition 50 C G E A N T t r a c k i n g code 56 D E x t e r n a l C o n v e r s i o n code 60 E Example Results 68 iv List of Figures 1 Feynman diagram for 7 emission 6 2 Feynman diagram for 7 exchange 6 3 Radial distribution of event vertices 20 4 Z distribution of event vertices 20 5 Radial distribution of event vertices; dotted line shows result of modelling detector acceptance 6 21 Z distribution of event vertices; dotted line shows result of modelling detector acceptance 7 21 Simplified G E A N T program structure (taken directly from the GEANT User's Manual) 8 31 Total cross section for pair production (barn/atom) as calculated using G E A N T parameterized fits of cross section data 9 40 Total cross section for Compton scattering (barn/atom) as calculated using G E A N T parameterized fits of cross section data 40 10 The S I N D R U M detector 51 11 Real data distribution of the "missing" E~, energy. The ne e~ peak is + on the left, the Dalitz peak on the right. The width of the ne e~ peak + reflects the resolution of the S I N D R U M detector. The fact that the peak is not centered at 0 is due to imperfect knowledge of the magnetic field. 70 12 x-distribution for real and simulated data 71 13 y-distribution for real and simulated data 71 v Acknowledgement I would like to express my graditude to the National Science and Engineering Research Council ( N S E R C ) for financial support, enabling me to pursue my studies; thanks also to my supervisor Dr. C. E . Waltham for his many hours of patient discussion and helpful suggestions, and to the members of the S I N D R U M group, most especially Dr. A . van der Schaaf, Dr. W . Bertl, and M r . C. Niebuhr, for their help and suggestions. vi Chapter 1 Introduction The collection of the theories of weak interactions, quantum electrodynamics ( Q E D ) and quantum chromodynamics ( Q C D ) known as the Standard Model has been very successful in explaining experimental results in particle physics over the last several years. The Standard Model is known to be far from complete, however, and in order to gain insight into its failings, physicists are continually testing its predictions. Suppressed in conventional theory, rare decays are sensitive to second-order effects and delicate, model-dependent cancellations, thus providing an opportunity to test higher order terms of the Standard Model. Rare decays also serve to provide limits for nonconventional theories by virtue of their sensitivity to "extra" physics. However, their very rareness makes them difficult to study experimentally. Measurements of rare decay parameters are high-precision experiments, requiring a complete understanding of both the apparatus involved and all possible backgrounds to the signal. Simulations are necessary both in the analysis and in the experimental design stages. W i t h a Monte Carlo simulation, one has the advantage of being able to examine "data" before and after it has passed through the detector, thus enabling a clear understanding of the action of the apparatus. During the design of the experiment, it is necessary to understand the signature of the decay in all areas of phase space, and to tailor the experimental apparatus so that it does not cut into the signal. For understanding of the background, all possible processes aping the decay signature must be understood, and the apparatus designed to cut them down as far as possible. Simulations also 1 Chapter 1. Introduction 2 provide an estimate of the systematic errors involved in the experiment, errors due to "imperfect" performance of the measuring apparatus. In the ideal case statistical errors are adjusted to be less than the systematic errors, and so the experimentalist can formulate a clear idea of how much data is required to extract the signal. In the final analysis, the simulation is used to provide a theoretical fit to the data. A l l background contributions to the signal are added together, together with the simulation of the decay of interest, and the result is fitted to the data, using the parameters to be measured as fit parameters of the Monte Carlo. The Monte Carlo simulation is an important tool in understanding both the experimental data and the detector, and its function can thus be divided into the two broad categories: 1. modelling the data, and 2. modelling the detector response. This document will outline the general procedures used in writing Monte Carlo code, and is divided into the above headings. Illustrations are provided from rare 7r° decays. Once each event has been generated, it is passed to the detector simulation part. Usually the tracking is left to software packages specially designed for the purpose of detector modelling, and this document will describe one such useful package called G E A N T , developed at C E R N . Appendices A through D list coded examples, and some final simulation results are shown in appendix E . Chapter 2 Data Simulation: Event Generation At the heart of any Monte Carlo simulation lies the event generation. It is here that the event kinematics are chosen and that the event is made 'real' by transforming it into the laboratory frame. Once the event has been generated, the relevant kinematical quantities are passed to the tracking routines and the detector response is simulated. From theoretical considerations a differential cross-section da(xi, x ,..., 2 x) n can be calculated, as a function of the kinematical quantities ,T,-. The cross-section expresses the probability of seeing the final state characterized by the x . % Once we have an expression for da, Monte Carlo sampling methods are used to choose random values of the kinematical quantities according to this probability distribution. Once the primary kinematics x± are fixed, the event is worked out in detail, any arbitrary orientation angles are chosen, and the appropriate Lorentz transformation is applied to bring the event into the lab frame. The event can then be passed to the tracking routines of the Monte Carlo, where the decay products are traced through the simulated detector, and its response is modelled. Secondary decays, showering, bremsstrahlung, pair production, and other such processes may be simulated at this point. The following sections describe the outlined general method of coding event generation in more detail, and will be illustrated by examples of 2-, 3-, and 4-body decay of the TT°. 3 Chapter 2. Data Simulation: Event Generation 2.1 4 T h e Differential B r a n c h i n g R a t i o In order to model any reaction we must be able to predict in a quantitative way the probability of all possible outcomes. In particle physics, this means that if we wish to simulate, for example, the reaction 7T~ 4- p —r n + e + + e~ we must be able to predict the energy distribution of the n, the opening angle distribution of the e e" pair, etc. This information is given in a quantity called the "differential + cross section" dcr(xi,X2, • • •, x ), 1 n where the X{ are the (independent) kinematical vari- ables (energies, opening angles) which describe the reaction. The "total cross section" obtained by integrating da(x) over all X{ is simply the probability for the reaction to occur at all (relative to some other process). We note that this means that: J •••J da(xi,..., x) n dxi • • • dx n — total cross section (2-1) to within a normalization factor. Usually the explicit calculation of da(x) is left to theorists who are more familiar with the intricate workings of the required models, but a brief outline of the general principles involved follows. More explicit details can be found in most quantum mechanics textbooks; for example, Cohen and Tannoudji [1]. The summary given here is taken principally from Halzen and Martin [2]. Using perturbation theory and the Schrodinger equation, we can find the eigenstates of a system which has been under the influence of a small, transient potential V. Initially, we have the system in a state <fo, acted upon by an unperturbed hamiltonian Ho for which we know the eigenstates (the <f> orthonormal): n Ho^n = E (f) n 1 n o r "differential branching ratio" in the case of a decaying system Chapter 2. Data Simulation: Event Generation 5 The final state is then given by a linear superposition of the eigenstates of the unperturbed system: To find the coefficients a we multiply through by (j>*j and integrate, using the orthonorn mality condition. The coupled differential equations thus obtained can be solved by making the approximation that the initial conditions hold true for all times. We obtain a (t) = -i J* dt' j d x <f>) V fa '( /-SO*' s 3 f e to first order, where the perturbation is turned on at time —T/2. Thus, when the perturbation is turned off at time T / 2 , we have T fi = a (T/2) f = -if dt [ J-T/2 = d x[4>j(x)e- ^}*V(x,t)[U^y- } 3 J iE lEtt -i J d x<t>){x)V{x)(j) {x) A l in covariant notation; which defines the "scattering matrix" T/,-. 2 We can extend the above technique to the realm of relativity, where the correct wave equation to use is the Dirac equation (using covariant notation) If we consider, for example, the scattering of an electron by an electromagnetic potential, we consider Maxwell's equations and make the substitution A< _» y< p 7 + A» e T h i s scattering matrix is related to the transition probability from state i to state / per unit time, Wfi, which figures in Fermi's Golden Rule, by 2 Wu = l i m Chapter 2. Data Simulation: Event Generation Figure 1: Feynman diagram for 7 emission 6 Figure 2: Feynman diagram for 7 exchange so that the Dirac equation becomes with the perturbation Carrying on as in the Schrodinger case, we obtain T = -iJ jfA^x fi (2.2) with defined as the "electromagnetic transition current" between the initial and final states. Note that the wavefunction here includes both momentum and spin, i.e., more ex- plicitly, j^-eufWie^'-rt* (2.3) where u is a Dirac spinor. Figure 1 shows the corresponding Feynman diagram. Now we can calculate the scattering of the same electron by another electron (or any other charged particle). We identify the source of the electromagnetic field A^ as the second electron, and solve for the field: • A " = ft 2 Chapter 2. Data Simulation: Event Generation 7 with solution Q where q = p — p , the momentum carried away by the photon. Figure 2 shows the a c corresponding Feynman diagram. Substituting the result above into equation 2.2 leads to T = -iJ*xjl (-^j fi to first order. Performing the x integration, (recall the explicit definition of the transition current j , equation 2.3) we get M H T = (-eM 7 u )(27r) 8(p + p - p - p ) (-euc-y^a) — -i (27r) 8(p 4 a D M 4 i a b c d +Pb - Pc~ Pd) M The delta function serves only to conserve 4-momentum, and M. is called the "invariant amplitude" or "matrix element", and is where the interesting physics resides. The 1/q 2 factor is known as the "(photon) propagator". The experimentally useful quantity, the cross section, is defined by no. final states cross section = da = Wu —. . initial flux and hence can be related to the matrix element: do = V^-dQ r (2.4) where F and dQ symbolically represent the incident flux and the Lorentz invariant phase space factor; 3 F= \v \2E -2E a a b where Pa v = —f- = velocity of beam particle a 3 for an excellent, more in-depth discussion of this important concept, see reference [12] Chapter 2. Data Simulation: Event Generation 8 and d Q = II ( 2 ^ : ( 2 7 r ) < ^ + W - P i - • • • - Pn) 4 (2-5) for a + b —> 1 + 2 + . . . + n. In the case of particle decays, we cannot define a cross section since there is no initial flux. We speak instead of a "differential decay rate" ("differential branching ratio", "partial branching ratio") dT where for a —> 1 + 2 + . . . + n. In summary, then, da (dT) is what we will want to sample. The physics is contained in M., which is given by adding up the relevant Feynman graphs in perturbation theory. The above discussion describes the basic premises of quantum electrodynamics (QED), which describes interactions between two charged point particles. In the case of hadrons, we are not dealing with point particles, but with extended structures. We can either reduce the mesons to their constituent quarks (which are charged point particles) and apply Q E D , or approximate the hadron as pointlike, modified by "form factors" which take into account the meson's internal structure. While the quark model looks to be the more elegant and fundamental, it is not any more complete than the form factor approach. Two complications arise: 1. the quarks have an extra internal degree of freedom, called "colour", so that their wavefunctions must be expressed in momentum, spin, and colour space. As mesons have no net colour (they are "colour singlets"), when describing them, we must sum over all quark colours. 2. the quarks inside the meson axe bound, not by electromagnetism, but by another force called the "strong force", for which a candidate theory is quantum Chapter 2. Data Simulation: Event Generation 9 chromodymamics (QCD). Since Q C D is not as calculable as Q E D , the usual technique is to introduce a form factor, where this time it serves to correct for the binding of the quarks. We quote here an example result of such a calculation (references [3], [4], [5]); dT for the process ir° —> e + e~ + 7. + - -c£*/;*^K -3 ( » ^) i r , i + , + (2 - 6) where m is the electron mass and / i is the mass of the 7r°, and 4?n 2 ~x^~ Here, f(x) is the form factor discussed, and, to first order in its x expansion, / is of the form f(x) = 1 + 2ax with a some constant. The equations derived are parametrized by x and y, where x — invariant mass of the e e~ pair y = 2.2 + energy partition = -rz, zri \P++P-\ (2-7) (2-8) Monte Carlo Sampling Methods The kinematics of each event are chosen randomly according to the differential branching ratio. This is done either by straightforward sampling methods or by using specially designed software. A n excellent review of the general problem of Monte Carlo sampling, as well as a discussion of Monte Carlo phase space integration, is given in [6] and [7]. The simplest way of generating random numbers according to a specific distribution is by inverting the probability distribution function and hence transforming random Chapter 2. Data Simulation: Event Generation 10 numbers on the interval [0,1) directly to the desired distribution. For example, to generate random numbers x distributed according to the exponential function e~ , we ax solve the integral equation u= fe- 'dx' Jo a x where u is a random number on the interval [0,1), for x, to obtain x = l n ( l — au) a Most functions, however, are not analytically invertible, so more indirect approaches are necessary. Given any distribution function F(x) one can choose random values of x which follow the distribution by the following rejection method. A random value a-'o is chosen, where x lies in the domain of f(x). F(x ) is evaluated, and then a new random number 0 Q y is chosen, where y lies in the range of F. The random number x is kept if y < 0 Note that this is an exact technique. F(x ). 0 It is adequate for functions which are slowly varying over their domains. For more sharply varying functions, the above technique made more efficient by dividing the domain of F into smaller regions. Each region Ax{ has associated with it a cumulative probability pi, j F(x) dx = pi with r^max I F(x) dx = l Note that this means that we must know the integral. The first random number generated decides i, the region to consider, by random number < p; Once the region has been decided, the preceding technique is applied. Chapter 2. Data Simulation: Event Generation 11 This technique can be extended to multi-dimensional integrals, even if they are not separable. The general idea is to first generate Xi, then use it to generate x , and so 2 on. We again need the cumulative probabilities, integrating over the quantities we are not choosing: Once X i is fixed, we calculate the cumulative probabilities for the next level down, and repeat. While the idea is straightforward, the programming is usually not, and numerical problems may arise when there are many dimensions, and the function is sharply varying. The problem becomes one of dividing up the integration space in order to optimize the performance of the numerical integration routines used to calculate the cumulative probabilities. The situation becomes even more complicated when parts of the integration space are cut out by introducing extra constraints on the limits; it is often not necessary to generate random numbers over the whole space, as detector acceptance cuts away events. For these reasons, one usually turns to packaged software. A n important software source is the large accelerator facility C E R N (Centre E u ropeenne pour la Recherche Nucleaire) at Geneva, which employs an entire division of computer specialists dedicated to developing and supporting programs of use for particle physicists. Programs and documentation in the C E R N Computing Center Program Library are made publicly available. Extensive useful scientific literature is 4 also available from C E R N ' s Scientific Information Service. DIVON C E R N has acquired and supports a program package called D I V O N which performs both integration and sampling. It is invoked by a straightforward calling sequence, further information obtainable from: Program Library, Division D D , C E R N , 1211 Geneva 23, Switzerland 4 Chapter 2. Data Simulation: Event Generation 12 see [8] for the user's guide and [9] for a detailed description of the algorithm. A brief summary of D I V O N ' s capabilities and of the partitioning algorithm follows. See appendix A for a coded example. D I V O N partitions, then integrates and/or samples a user-defined function F, which can be multidimensional. Absolute (numerical) limits of integration are passed to D I V O N by the user. Further interdependent limits and cuts may be coded within the function. Various common blocks allow for user control of accuracy and termination of the partitioning and integration. 1. P a r t i t i o n i n g : D I V O N divides the integration region into hyper-rectangular subregions Ri, each characterized by its variation V^-, the range of extreme values of F over the region. The extreme values are found by a numerical maximization/ minimization procedure, which uses as a starting value the maximum/minimum function value of N P O I N T 5 random, uniformly-distributed vectors drawn from Ri. For proper partitioning, it is recommended that N P O I N T be put to a high value, say, 150-200. The spread Si = \V{ -vol(Ri) is a measure of the Monte Carlo integration error within Ri and hence gives an estimate of the contribution of Ri to the total integral uncertainty. The partitioning algorithm selects the region with the greatest spread and subdivides it, repeating until /n \ 1/2 S = \ Y, i) S ^ SPRDMAX where S P R D M A X is user-defined, or until the number of function evaluations exceeds the user-specified limit M A X P T S , whichever comes first. S P R D M A X should be in the range of 0.2-0.5. Depending on the difficulty of the integration, this will require on the order of 10,000-100,000 function evaluations. In the case of Dalitz decay 7r° —> e + e~ + 7, integrating equation 2.6 requires 10 function + 5 COMMON/SAMPLE/NPOINT s Chapter 2. Data Simulation: Event Generation 13 evaluations. After partitioning, the variation of F over each subregion is small, allowing for accurate numerical integration. 2. I n t e g r a t i o n : D I V O N supports several alternative numerical integration procedures: • quasi-uniform Monte Carlo estimate (recommended) • pseudo-random estimate with importance sampling • 2 , 3 , or 5 nd rd t h degree quadrature estimates • a user-provided algorithm. Accuracy of integration is controlled by the variable N P T , the number of samples drawn in each of the partitioned regions. The user should experiment with the value of N P T to obtain the desired accuracy. N P T values of 10,000 are not unusual. 3. S a m p l i n g : D I V O N keeps an integral estimate of each subregion as it partitions, and uses this as a measure of the relative probability when choosing random vectors. In order to increase the accuracy of these integral estimates, it is recommended to force D I V O N to integrate the entire region before sampling. Two sampling options are available: • G E N P T ( N , X , W T ) , which fills N-dimensional random vector X , weighted by W T . This results in a distribution statistically equivalent to a sample distributed according to F. • R A N G E N ( N , X ) , which fills N-dimensional vector X exactly according to the coded distribution function. Chapter 2. Data Simulation: Event Generation 14 In order to test the results obtained with D I V O N , a good qualitative idea of the behaviour of the integrand is necessary. The first thing to check should be the partitioning; a useful technique is to dump the results of the partitioning using the facility B U K D M P , and then to use whatever graphics routines are available to plot them. It is then easy to see if D I V O N is producing a reasonable partitioning—the most sharply varying regions of the function should be heavily partitioned. Note that D I V O N is not designed to recognize or exploit any symmetries inherent in the problem, so that it may be necessary to adjust the integration limits, to divide the problem up into several regions, in order to achieve the desired results. The total integral should be checked against any theoretical (and/or experimental) total cross sections (c.f. equation 2.1). A quantitative knowledge of the integrand and integration or of closely related functions is also invaluable in testing. It is recommended to test the integration in a cumulative way, in order to test the results over the entire domain. P h a s e Space Integration It must be noted that some theorists calculate only the model-dependent and hence more interesting matrix element, leaving the phase space factor dQ (equation 2.5) to be handled by those wishing to use their results in a Monte Carlo simulation. In other cases, theorists have not been able to provide a matrix element calculation (in the case of decays forbidden in the Standard Model, where the details of the interaction are not known, for instance), so that it may be necessary to model the process using a constant matrix element, or one which is very simple in nature. In nuclear physics, also, experimentalists frequently use phase space Monte Carlo calculations modified (by a weight, dependent on the 4-momenta of the particles involved) to simulate resonance formation or a momentum-transfer-dependent interaction. Many phase space sampling programs exist, and the literature is extensive. The Chapter 2. Data Simulation: Event Generation 15 programs are essentially Monte Carlo integration routines dedicated to integrating and sampling dQ for n up to about 20. R A M B O and F O W L , ([10], [11]) both available from the C E R N program libraries, are two such packages. F O W L is a stand-alone program and therefore more useful as a "first-estimate" (non-detector-specific) Monte Carlo than as an event generator. R A M B O is organized as a subroutine call and hence more amenable to use as an event generator, but it does not have the options or flexibility of F O W L . It is of course possible to use D I V O N in this context, simply by supplying it with the required phase space integral. Reference [13] gives suitably parametrized integrals for 2-5 body phase space. 2.3 Laboratory Frame Kinematics The random kinematical quantities having been chosen according to the differential branching ratio, the event can now be materialized in the laboratory. From the chosen quantities, the entire event must first be reconstructed in the frame in which the branching ratio calculation was done. For particle decays, this is naturally the rest frame of the decaying particle, while for particle collisions, it will be the centre-of-mass system. It is possible that neither system corresponds to the laboratory frame, and so the event must be Lorentz-boosted into the lab. For 2-body reactions, knowledge of the initial conditions fixes the event kinematics but for two arbitrary orientation angles, while for 3- and many-body reactions, there will be freedom to choose opening angles and particle momenta. The above generalities become clearer once applied to a specific system. Consider the following situation: a beam of 7r~'s is incident on a target containing only protons (liquid hydrogen). The 7r~'s are slowed down by passing them through a lead sheet ('degraded' or 'moderated') until they are virtually at rest when they interact with the Chapter 2. Data Simulation: Event Generation 16 protons, also at rest. We are interested in resulting e e pairs. We have the reactions: + ir~+p -> 7r° + n (2.9) with 7T —> e + e (2.10) -»• e + e~ + 7 (2-11) + e + e~ + e + e -» + + _ (2.12) and also 7r-+p -» e + + e~+n (2.13) producing various e e~ combinations. Consider reaction 2.9 above. Knowing the initial + conditions, it is straightforward to derive (we neglect the binding energy of the ir~p system and any kinetic energy in the initial state): E„o = 137.8556 M e V fa = 0.2037 E = 939.9724 M e V = 0.0299 n (3 n The kinematics is now completely fixed and we need only orient the event in space. We do this by choosing the direction of either the 7 r ° or the neutron. This is done by choosing p , x p y and p randomly, which, expressed mathematically, means sampling z the distribution dp dp dp x subject to the constraint y z Chapter 2. Data Simulation: Event Generation 17 Converting to polar coordinates, this is equivalent to sampling J j d<f>rf(cos6) which corresponds to choosing (f> randomly from 0 to ISO degrees and cos# from -1 to 1. Reaction 2.10 is handled in a very similar manner. The kinematics is first fixed in the rest frame of the 7r°, and then the entire event is Lorentz-boosted into the lab frame: E lab = jE -j/3 -p Pub = p-MiE-ft-p- 1 7+ 1 1 with P PIT = 0 Eir° 1 7 the usual Lorentz quantities. For reaction 2.13 sampling is performed to obtain the kinematics from the differential branching ratio. In this case, we use the results from Kroll and Wada [3] who express their differential branching ratio in terms of x and y, see equations 2.7 and 2.8. Once these have been chosen, we may determine the event kinematics: El -m + x 2 t 2 n ^pair — Q jp 2E t to E+ — ~^(Epair ~\~ Impair ' J/) E- = Ep i — E+ a r Chapter 2. Data Simulation: Event Generation with the angle between the e 18 and e given by + 2E,E_ - 2m - x 2 2 e cos a — where E ot = E - + E = m - + m t T v p p (again neglecting binding energy). We now align the neutron along the z-axis and "glue" the e e~ pair to it, allowing + one random angle in the x — y plane; we have the freedom to spin the e e~ fork about + the neutron direction. Once we fix this, we pick a random direction for the neutron and realign the entire structure by the 8 and <f> angles of the neutron: Pe,new — Pe,old where TZ is the rotation matrix cos 9 cos <f> n TZ — — sin 0 n n cos 9 sin <f> cos 8 n y n — sin 6 n n 0 sin 6 cos <j> \ n n sin 9 sin <f> n cos 9 n n j Reactions 2.11 and 2.12 are modelled as above, and once the kinematics are fixed in the 7r° frame, the events are boosted into the lab. Testing the kinematics is straightforward; one tests the final laboratory event for momentum and energy conservation. 2.4 The Event Vertex To fix the event completely we must localize it inside the detector. This involves generating the event vertex, the coordinates where the interaction has supposedly taken place. Usually this will be inside the target, if we are modelling a fixed-target experiment. It is necessary to model a reasonable distribution of event vertices; not all events Chapter 2. Data Simulation: Event Generation 19 will be generated at exactly the same location. The effects of this vertex distribution may not be negligible; the resulting particles will have different amounts of target material to pass through and hence experience different energy losses. It is therefore usually desirable to model this stop distribution rather accurately. In most cases, one postulates a particular distribution (for example a box or gaussian) and uses this as a first approximation. If more accuracy is desired, one looks at the real data and either fits or digitizes the actual vertex distribution, and uses the resulting parameterization in the simulation. Digitization is done by simply taking a histogram of the stop distribution and using its (normalized) bin contents to calculate the cumulative probabilities which will be used sampling the event vertex. Note that having the input of the simulation the same as the output of the detector does not guarantee that the output of the simulation will match the data! Modelling the stop distribution is thus an iterative procedure; to account for the effects of detector acceptance on the measured distribution, one must keep making adjustments to the parameterization until the output of the simulation matches the data. Figures 3 to 6 show examples of how the procedure works. In figures 3 and 4 we show the real vertex distribution taken from experimental data, on which we superimpose the results of the digitization. The digitized distributions are used to generate events, and, after analyzing these simulated data, we plot the resulting event vertex once more. Hence the difference in the distributions shown in figures 5 and 6 are due to the effects of detector acceptance. For most applications, the distributions shown are accurate enough. 2.5 M o d e l l i n g E x t e r n a l Conversions If the final state of the reaction of interest involves electrons, then external conversions can be an important source of background. External conversions occur when photons er 2. Data Simulation: Event Generation solid line: real doited data line: d i g i t i z a t i o n result vertex R (mm) Figure 3: Radial distribution of event vertices I I i solid line, dotted real data line: digitization result -120 Figure 4: Z distribution of event vertices Chapter 2. Data Simulation: Event Generation 21 .025- solid line: real data dotted line: monte carlo .020 • 2f> vertex R (mm) Figure 5: Radial distribution of event vertices; dotted line shows result of modelling detector acceptance. solid line: real dotted data line: m o n t e carlo T3 <U £ .016 10 E i— o —> .010 c o 3 vertex Z (mm) Figure 6: Z distribution of event vertices; dotted line shows result of modelling detector acceptance. Chapter 2. Data Simulation: Event Generation 22 interact with matter, resulting in the production of electron-positron pairs. Several mechanisms exist for electron production, each with very different characteristic cross sectional behaviour, which will determine its importance as a source of background. Usually the dominant processes are direct pair production 7 —> e + e~ and Compton + scattering (the scattering of a photon by an atomic electron, which can be knocked free in the process). Compton scattering is "forward peaked", i.e. the electron has a high probability of carrying off most of the photon's momentum and thus ends up travelling in the same direction as the incident 7. Electrons generated by either of these processes are of course unrelated to the reaction being studied, but may pass through the trigger system and be mistaken for the real signal. They must therefore either be removed from the data by improving the selection criteria, or be simulated and included in the final fit. In either case, modelling is called for. In order to cut down on external conversions, detectors are usually made as lowmass as possible. The probability for a photon converting is then very small, making the process extremely time-consuming to simulate in the direct manner. Nature, of course, has no problems with this, being able to generate many more photons than our computer can generate random numbers. Thus, instead of following each 7 through the detector and waiting to see if it will undergo conversion, the technique is to force 6 it to convert at a given location, and then to calculate the probability of the conversion having occurred at that point. This probability is then stored along with the usual event kinematics. The location of conversion is varied. In the final analysis, the probability is used as a weight for the event; in other words, instead of simply counting the number of events with certain characteristics and plotting this count (histograms), one increments histograms weighted by the probability. i.e. stepping the 7 through the detector, at each step calculating the probability for conversion p and then allowing conversion if (random number) < p . Obviously if p is small we will need many 7's to see an event. 6 c c c Chapter 2. Data Simulation: Event Generation 23 Other techniques are sometimes employed to model these rare processes. One can simply increase the probability for the interaction to occur. This is dangerous, however, because second-order effects are then affected in a non-linear way. For example, if external conversions go as e, then the probability for a photon to undergo two conversions goes as e . If one now increases e by a factor of 10 to make conversions more likely, 2 the probability of second-order conversions is increased 100-fold, which is non-physical. Another technique used is to follow the steps outlined, forcing the conversion/scattering at a particular location, but then instead of sampling the appropriate cross section to find the momenta of the final state particles, one forces a particular final state and carries along its probability as yet another weight. This is useful for distributions with long tails, where only a small part of the phase-space is interesting. If the direct sampling approach is used in such cases, the computer will spend most of its time in the tails. The technique of forcing an event and subsequently carrying along its weight can be used for any process which has low probability. In the case of 7r° decay, for instance, the decay modes involving electrons are rare: 7T° - f e + e~ + TT° - > e + e~ + e + e~ (1.198 x 10" branch) + + 2 7 + (3.24 x 1 0 - 5 branch) (2.14) (2.15) Instead of generating a horde of 7r°'s and waiting for one of reactions 2.14 or 2.15 to happen, we force each of these reactions and add them up at the end, remembering that 2.14 contributes about 185 times more electron pairs than 2.15. Since a more detailed discussion of the techniques of modelling external conversion and Compton scattering requires knowledge of particle tracking, we defer further discussion of these processes until after the introduction to G E A N T . Chapter 3 Detector Simulation: Tracking Once the event kinematics have been fixed according to theory, the individual particles must be tracked through the detector and its response simulated. The simulation must perform the following: • S t e p p i n g : each particle produced by the theory and kinematics subroutines must be stepped through the detector. This may involve simulating the particle tracks through a magnetic field. It will involve stepping the particle through matter, in which case the probability of it undergoing some kind of interaction must be taken into account. Scattering, showering, annihilation, secondary decays; all of these and more may be possible and, since they affect the 4-momentum of the parent particle, they will affect detector response (see below). Such processes may also produce daughter particles, which can then cause spurious signals. • Digitization: as well as insensitive matter, the detector contains sensitive elements (wire chambers, Cerenkov detectors, calorimeters, etc.) which produce signals when particles traverse them. Their response depends on the nature of the particle; its energy, its charge, its direction of flight. For correct simulation, the location and behaviour of the sensitive elements must be known. This is a complex task, for in real life, the resolution of the detector is a function of its design (spacing of wires in drift chambers, scintillator thickness, etc.) and the 24 Chapter 3. Detector Simulation: Tracking 25 intrinsic resolution of the electronics. The two cannot be separated when examining data from the experiment. The design resolution can be understood in terms of the geometry and materials used, but to model the response of every piece of electronics is a time-consuming taks. For this reason, the usual approach when writing the simulation is to assume perfect electronics, and then to smear the final results until the data is matched. If more accuracy is desired, the noise of the electronics can be more accurately simulated. • O u t p u t : ideally, the simulation should produce data identical in format to the real experiment. The output of the simulation must look indistinguishable from what comes out of the detector's electronics. Then both the Monte Carlo data and the experimental data can be analyzed in exactly the same manner. • G r a p h i c s : while not necessary, it is exceedingly useful to be able to see directly what has been programmed into the simulation. It is most useful if the simulation provides both views of the detector, so that mistakes in positioning the various elements are easily spotted, and views of the particle tracks, which is of help when debugging the tracking routines. The would-be author of any simulation must hence be familiar with the detector to be modelled and the principles of operation of its various sensitive elements, as well as knowing the specific format of the output of the electronics and a clear idea of how the final experimental data will be written on tape. For small detectors with few parts and a straightforward geometry, no magnetic field, and no calorimeters (which require showering code), it is not too difficult to write a simulation fulfilling the above criteria. For anything more complex, however, the programming quickly becomes very involved. In answer to the increasing demand by the large experiments clone at C E R N for such simulations, their Data Handling 26 Chapter 3. Detector Simulation: Tracking Division has developed an all-purpose detector simulation package called G E A N T . The following sections will present a general introduction to G E A N T , but are not intended as a substitute for the very thorough and excellent G E A N T User's Guide [14]. We will attempt only to give the reader a general idea of the use of G E A N T , and will mention known bugs (which the User's Guide does not) and stress some of the aspects of G E A N T we have found to be of importance. Coded examples can be found in appendices B and C. 3.1 G E A N T : Introduction and General Philosophy G E A N T is based on the memory-management system Z E B R A (also developed by C E R N ) which organizes information into arrays or "banks", which are then addressed by pointers. The pointers are carried in G E A N T common blocks, which permit the various subroutines to communicate with one another when manipulating the stored data. G E A N T routines are basically an interface between Z E B R A and the user, with a specific application in mind; the detector definition routines allow the user to enter information into the Z E B R A banks, the tracking and digitization routines manipulate the stored information. Thus, the user need not be familiar with the workings of Z E B R A when using G E A N T . G E A N T comes equipped with a graphics facility, which is based on C E R N ' s G D 3 / T V software. Interfaces to several other graphics packages are also available in principle, using special linking sequences; if G E A N T is available on your computer system your system manager should be able to provide these details. G E A N T is designed as a series of subroutine calls within a specific framework. Most of the subroutines are general-purpose, automatically invoked, and hence "invisible", but G E A N T allows the user to control the tracking and digitization at key points by calling user-defined subroutines. In general, G E A N T subroutines begin with the letter Chapter 3. Detector Simulation: Tracking 27 ' G \ while the user-written routines called by G E A N T start with ' G U ' or ' U G ' . We present here a general program outline, and in the sections following give more details on the various user-written parts. It is not necessary to supply G E A N T with every user routine; G E A N T is equipped with default user routines which perform the minimum actions required to run the program. The reader is referred to the G E A N T manual for full details on calling syntax, as well as a complete description of the action of the G E A N T subroutines and of the default values of all parameters. We also do not describe here the graphics capabilities of G E A N T . Note that G E A N T works in the following system of units: centimetre, second, kilogauss, GcV, GeV/c, degree. The main driving program of G E A N T must be provided by the user. This main program must allocate the memory space for Z E B R A by a call to G Z E B R A . Further, M A I N controls the three phases of G E A N T : • initialization; in the form of a call to U G I N I T , • event generation and tracking; by a call to G R U N , and • termination; by calling U G L A S T . Each of these "phase control" subroutines, whether user-written or not, calls a number of subsidiary routines which perform the actual data manipulation. A more detailed summary follows. 1. I n i t i a l i z a t i o n : U G I N I T . This user-provided subroutine sets up various G E A N T structures which are necessary for running. This involves a number of subroutine calls (some of which are optional): • GZINIT—allocates memory for the dynamic structures which will contain the detector and tracking information (tailors Z E B R A for G E A N T use). Chapter 3. Detector Simulation: 28 Tracking © GINIT—initializes G E A N T common blocks with default options (including switches for various physical processes, the initial value for the random number seed, the number of events to generate). 1 e GFFGO—reads any desired input data cards to override the default options set in GINIT (a facility exists for the creation of new options and associated data cards). A complete list of available data cards is given in the manual. e GDINIT—opens the drawing banks for graphics (if desired). • TVBEGIN—initializes the graphics package (if desired). G E A N T also needs to initialize its tables of particle properties, material properties, and cross section tables for modelling physical processes: • GPART—loads particle lifetimes, charges, masses, etc. • GPHYSI—calculates interaction cross sections, energy loss tables, etc. • GMATE—sets up properties of the most commonly used materials (radiation lengths, density, atomic number). We call the user's attention to the fact that the numbers given in the G E A N T manual for the materials hydrogen, deuterium, helium and neon are for the liquid state. The appropriate values for gaseous helium.are: density = 1.78 x 10" g cm' radiation length = 1.075120 x 10 cm absorption length = 6.80 x 10 cm 4 3 6 5 N o t e that while the G E A N . T manual states that the default switch for bremsstrahlung is I B R E M = 2 , that is to say, bremsstrahlung is by default calculated but the resulting gamma ray is not tracked, this is not the case. The gamma is in fact tracked ( I B R E M = 1 ) , and if this is not desired, the user should include the data card B R E M 2 in the input. 1 Chapter 3. Detector Simulation: 29 Tracking G E A N T uses these material properties to determine the step length when tracking particles through these media. We mention here that a bug exists in G E A N T ' s use of the vacuum. G E A N T calculates the step length using 1/density, which will give a very small number for vacuum rather than a large step as is desired. The user should therefore overwrite the density of vacuum with a small number, say, 1 x 10~ . 6 • GSTMED—defines the tracking media. A tracking medium is defined by, among others, its material (previously defined in G M A T E ) , a flag for the presence of a magnetic field, and parameters specifying the tracking precision. Subroutines for overwriting or adding to the above tables exist, and, if used, should be called here. The user may also define mixtures and compounds of previously defined materials at this point. Finally, as a last step in the initialization, the user must define the detector which G E A N T is to model: • user code—sets up the detector by specifying materials used and detector geometry, also specifies which parts of the detector are sensitive. We discuss the techniques involved in section 3.1.1. 2. G e n e r a t i o n a n d T r a c k i n g : G R U N . This phase controls the processing of events. It is here that the (invisible) loop over events is located. The user selects the number of events to be generated by the input card T R I G number-of-events. G R U N will call the following G E A N T subroutines: • GTRIGI—initializes event processing • GTRIG—generates and tracks one event. During this phase, the user is handed control at several points: 30 Chapter 3. Detector Simulation: Tracking — GUKINE—user code to generate event kinematics. — G U T R A K — c a l l e d at the beginning of every track, allowing the user to affect the tracking of different particles in different ways. — GUSTEP—gives the user control after every tracking step. — GUSWIM—tracks the particle along one step in the magnetic field. The user can specify a simple one-component field, or can define an inhomogeneous field in the user routine G U F L D . — GUOUT—called at the end of each track, allowing the user to accumulate any desired output. — GUDIGI—models detector response for the entire event. A more complete discussion of the tracking procedure is given in section 3.1.2, and we discuss further the modelling of detector response in section 3.1.3. • GTRIGC—empty data structures after every event. 3. T e r m i n a t i o n : U G L A S T . The termination of G E A N T usually simply involves closing the data structures down by a call to G L A S T , and ending any graphics by calling T V E N D . Figure 7 shows a program simplified flowchart. 3.1.1 Detector Definition Here we describe, in a general way, the use of G E A N T for defining the detector setup. A coded example can be found in appendix B . This code must be called by U G I N I T during the initialization phase of G E A N T , as described in the previous section. The detector is defined by a series of volumes, defines an initial mother volume, positioned in one another. The user and positions daughter volumes inside this. These Chapter 3. Detector Simulation: Tracking 31 M A I N (user) — GZEBRA initialisation of ZEBRA ivstcm, dynamic core allocation — UGINIT (user) — — — — — — GINTT GFFGO GZINIT GPART GMATE 'user code' initialisation of GEANT3 variables interpretation of data cards initialisation of ZEBRA core divisions and link areas creation of the 'particle' data structure JPART creation of the 'material' data structure JMATE description of the geometrical setup, of the sensitive detectors creation of data structures J V O L U M , JTMED, JROTM, JSETS — GPHYSI preparation of cross-section and energy loss tables for all used materials — G R U N loop over events GTRIGI initialisation for event processing GTR1G event processing - GUKINE (user) generation (or input) of event initial kinematics - GUT REV (user) \— GTREVE (loop over tracks, including any secondaries generated) I - GUTRAK (user) l—GTRACK control tracking of current track |— GMED1A find current volume/tracking medium |— GTVOL loop over successive media seen by the particle — f GTGAMA/GTELEC/.... tracking of particle according to type r-GUSTEP (user) recording of hits in data structure JHITS and of space points in data structure J X Y Z - GUDIGI computation of digitisations and recording in data structure JD1GI - GUOUT output of current event — GTRIGC clearing of memory for next event UGLAST (user) I— GLAST standard GEANT3 termination. stop Figure 7: Simplified G E A N T program structure (taken directly from the GEANT User's Manual) Chapter 3. Detector Simulation: 32 Tracking daughter volumes can themselves be made into mothers by postioning further daughter volumes inside them. Note that G E A N T does not think further than two generations; operations done on "a volume and its contents" affect only the mother volume and its daughters, not its grand-daughters. Each mother and daughter volume has associated with it its own reference system. The tracking is done in the so-called Master Reference System, and is denned by the way the user represents the kinematics of the events (in G U K I N E , see section 3.1.2), which will usually be some orientation of a cartesian co-ordinate system. This Master Reference System coincides with the local reference system of the initial mother. G E A N T provides subroutines for transforming from one reference frame to another. Each volume is specified by several parameters, which are passed to G E A N T when the volume is defined: • its S H A P E , one of several basic geometrical shapes available including cones, cylinders, boxes, trapezoids, polygons, spheres; • its dimensions, the number of which depends on the shape; • its number, which is a counter incremented in the order of volume definition; • its local reference system, which is defined by the shape used; and • its physical properties, depending on the tracking medium (as defined in the initialization) which fills it. G E A N T allows for both the positioning of new daughter volumes (with completely new parameters) inside a mother volume, and for the division identical volumes with the same parameters as the mother. of the mother into several Note that this division or slicing is allowed along only one of the axes of the local reference system of the mother. Chapter 3. Detector Simulation: Tracking 33 Local reference frames are defined for all the divisions, in such a way that the local origin is at the geometrical center of the new volume. Because of the proliferation of reference frames, the user should be well-armed with clearly dimensioned drawings of the detector to be simulated. It is at this stage that the graphics capabilities of G E A N T are most useful. 2 Volumes are denned using the routine G S V O L U , positioned using G S P O S , and divided using G S D V N . See the G E A N T manual for an excellent summary of the geometry tools. During the tracking phase, G E A N T steps the particle through each volume. When a boundary is reached, G E A N T must decide which is the next volume. Its first guess is then the volume with the next consecutive volume number. In the interests of faster tracking, therefore, it is important to define volumes in some logical order. It is also recommended that use be made of one of the optimization routines G S N E X T , G S O R D , or G S S E A R / G O S E A R . G S N E X T and G S O R D are called once after all volumes have been defined. G S N E X T inputs a user-defined list of volume numbers which specifies the order in which to search volumes during tracking (a default facility also exists, which limits the search to volumes with consecutive numbers). G S O R D causes G E A N T to create its own list for binary search during tracking by forcing it to check the ordering of volumes along a specified axis. For most applications, this routine is the most useful. G O S E A R should be called after every few events, and performs statistical updating of the search list 3 (only volumes "set" with G S S E A R are included in the list), so that G E A N T will test first those volumes which are entered most often. The final step in the initialization is to define sensitive areas of the detector for the W e draw the user's attention to the fact that G E A N T will not draw half-spheres, drawing full spheres instead. Special attention is thus required to ensure that these volumes are positioned correctly. G O S E A R outputs this updated list, but its output is rather cryptic; what is printed is the volume number and the number of times it has been entered, for each volume set by G S S E A R . 2 3 Chapter 3. Detector Simulation: Tracking 34 digitization phase of G E A N T . This is done using the routine G S D E T (and its associated routines). The digitization of the hits is a very experiment-dependent procedure, essentially consisting of modelling the electronics of the experiment. The user may therefore find it more convienient not to use G E A N T for the digitization phase, but to use home-grown routines to store and process the hits. In this case, it is not necessary to use G S D E T at all. It is possible to use G E A N T for temporary hit storage, and then to retrieve this information when digitizing; in this case, G S D E T is necessary in order to be able to access the hit banks later. 3.1.2 Tracking After initialization and definition of the geometry in U G I N I T , the M A I N program calls G R U N and enters the loop over events. Each event begins with its generation, by calling G U K I N E (see figure 7). In G U K I N E , the event kinematics are fixed, using the techniques discussed in the section on data simulation. The only inputs G E A N T requires are the momenta of all the ("primary") particles to be tracked, and the event vertex ( G E A N T supports more than one event vertex). This information is passed into the G E A N T data structures by using the subroutines G S V E R T and G S K I N E . Once the event has been generated, G E A N T passes control to the user routine G U T R E V to control the tracking at the event level. Here, for instance, counters can be initialized, user written event banks can be cleared, histograms relevant to whole events accumulated, etc. If nothing else, G U T R E V must call the G E A N T eventtracking routine G T R E V E . G T R E V E loops over all the particles it finds in the G E A N T data banks; this includes the primary particles passed by G U K I N E , as well as any secondaries generated during the tracking. G T R E V E will eventually call G U T R A K , allowing the user to take action at the beginning of every track. G U T R A K must call G T R A C K , which gets the particle type to be tracked from the G E A N T data banks, Chapter 3. Detector Simulation: Tracking 35 and calls G T V O L to step through each encountered volume until the track reaches its end ( i . e . leaves the initial mother volume). G T V O L accesses the data banks and finds the parameters of each volume encountered. The actual track stepping is clone by the routines G T G A M A , G T E L E C , G T H A D R , G T N E U T , G T M U O N , and G T N I N O (for each separate particle type), called by G T V O L . Here the tracking step size is calculated, depending on the present tracking medium, and the track is extrapolated by one step. If the current medium is in a magnetic field, G U S W I M is called, which will provide G E A N T with the field parameters required for the propagation of the track. G U S W I M is a user routine; its default action is to call • G H E L X 3 for one-component fields, • G H E L I X for fields tilted with respect to the reference frame, or, for inhomogeneous fields, • G U F L D , a routine provided by the user to find the components of the magnetic field at any point. The tracking-stepping routines also model any desired physical processes. The interaction lengths in the medium for various processes (pair production, Compton scattering, bremsstrahlung, in-flight decay, etc. See the G E A N T manual for a complete list of available processes) are calculated, and using these, the distances to the next interaction point are found. If this distance is less than the current step size, the track-stepping routine will call the appropriate interaction generator and create the interaction. The interaction generating subroutines sample relevant differential cross sections and produce "secondary" particles with appropriate momenta. 36 Chapter 3. Detector Simulation: Tracking The user routine G U S T E P is called at the end of every tracking step. Here the user must decide whether or not to keep the secondary particles generated in the last tracking step, by entering them into the banks by a call to G S K I N G . Once entered, they become part of the event and are tracked, until they pass out of the initial mother volume. In G U S T E P , the user can also abort tracking, if, for instance, the number of tracking steps becomes too large (the particle spirals endlessly in the magnetic field). In this routine, the user should include output options, in order to facilitate debugging of the tracking. The G E A N T routine G S C X Y Z stores the current tracking point in the banks, and using G P C X Y Z the user can then print out the tracking steps. Appendix C shows some example code for subroutines G U T R E V and G U S T E P . 3.1.3 Digitization and Output G E A N T provides the user with a set of data banks for recording "hits", the (userdefined) information regarding the track parameters (momentum, energy loss, coordinates, etc.) upon interaction with a sensitive volume (necessary for modelling detector response later), and for recording the final digitized data. G E A N T can not perform any digitizations itself; this phase of the simulation is essentially dependent on the types of detecting elements, electronics and on/offline computers used in the experiment, and will vary greatly from application to application. Hence G E A N T provides only the storage banks and routines to provide access to them, some general-purpose digitization programs for drift and multi-wire proportional chambers, as well as some routines to calculate track intersections with cylinders and planes. More detector- specific routines have to be provided by the user. If the user decided to use G E A N T ' s hit banks for storage, the routine G S A H I T can be used to store any desired hit information. G S A H I T should be called after every tracking step, i.e. in G U S T E P . G S A H I T can be called only for sensitive detector Chapter 3. Detector Simulation: Tracking 37 elements declared during initialization using G S D E T . If the user does not use the G E A N T hit storage facilities, provisions must be made in G U S T E P to determine and store any desired information regarding the sensitive volumes which have been hit. After each track has been followed to its end, G E A N T passes control to the user routine G U D I G I , where the detector response is to be modelled. Whatever structures have been used to store hits, they must now be unloaded and the information used to produce simulated experimental data. This will be clone either using G E A N T or by some user-provided routines. G E A N T provides the routine G F H I T S to unpack the information stored in the hit bank using G S A H I T . The actual digitization involves the translation of the stored hit information into simulated data; energy losses, track momenta, hit co-ordinates are transformed into drift chamber, scintillator, calorimeter signals. Other detector elements may be much more complex. Translating the hit energies and co-ordinates of drift and wire chambers, calorimeters, Cerenkov counters, etc. into the appropriate quantities obviously requires knowledge of how these elements function. The interested reader is referred to the extensive literature available as a full discussion of the operational principles of the detectors mentioned is out of the scope of this document. We consider here a only a simple hodoscope consisting of 64 scintillator strips. In real life, a "signal" occurs when an electron crosses the scintillator and, in slowing, loses energy in the form of ionization. The scintillator converts this ionization energy into photons, which are then carried through the plastic lightguides to the photomultiplier tubes attached to either end of the strips, resulting in a voltage pulse from both ends, proportional to the energy deposited. In the experiment, we will want to be able to record which of the 64 scintillators fired as well as the ouput voltages. The simulation must then also provide these numbers. Using G E A N T , we can define the scintillators as sensitive elements, and give them numbers. When one is traversed by an electron during tracking, we store 38 Chapter 3. Detector Simulation: Tracking its number and the energy deposited in it after every step. In G U D I G I we can then sum over all energy deposited and convert to the electronics' standard to simulate the signal. In the case of a scintillator, all that is needed is a "yes" or "no", depending on whether the energy deposited is above the scintillator threshold. Noise in the electronics can also be simulated at this stage. G E A N T provides a bank structure for storing the digitized hit data. The routines GSDIGI and G F D I G I store and fetch, respectively, from this bank. However, unless the user plans to use similar data structures for the experimental data, or to have different offline analysis programs for simulated and real data, these will not be too useful. It is recommended that the products of the simulation be as similar in format to the experimental data as possible, so that the same analysis programs can be run on both. This makes the task of understanding the final analysis program a much easier one; it also greatly simplifies its optimization and debugging. 3.2 E x t e r n a l C o n v e r s i o n s and C o m p t o n S c a t t e r i n g Now that the reader is familiar with the techniques used in stepping particles through the detector, we take up again the problem of modelling external conversions and Compton scattering. Before simulation of external conversion can be done, it must first be established where in the detector these reactions are most likely to take place and where they are to be simulated. Usually this will be in the most massive detector parts, any metal parts are good candidates: the target itself, the target housing/support, walls dividing the drift/wire chambers, etc. Locations that would result in obviously "wrong" event topology or outside of the sensitive detector area are not bothered with. Then, the 4 except in the case of rate calculations, where one is estimating the total number of events, regardless of topology, that the detector will see. 4 Chapter 3. Detector Simulation: 39 Tracking relative probability p for each of these sections to initiate a photon conversion must be x c found. This involves knowledge of the material properties and is of course dependent on the energy of the incoming photon. This information is contained in a quantity known as the "interaction length" A, the mean distance between two successive interactions. Extensive tables of material properties are given in [15] and the associated summary booklet, which are available free of charge by writing to the C E R N Information Service. As an example, we calculate the probability of external conversion in 2 cm of liquid hydrogen and compare this to the probability of conversion in 10 microns of mylar (C5.H4O2) f°r a photon of energy 50 M e V . In general, the probability p,- of a particle to undergo an interaction i is given by L Pi = T T = L- A — if the particle travels a distance L in a material characterized by density p, atomic number Z, and atomic weight A. NA is Avogadro's number, and <r,- is the total cross section for the interaction i to occur, and depends on the nature' and energy of the traversing particle and on the material. Figures 8 and 9 show plots of the total cross section for the process 7 —» e + e~ and Compton scattering in both liquid hydrogen and + mylar. Applying the formula above, we find that for 50 M e V photons, the probability of conversion in hydrogen is about 10 times greater than in the nrylar. Once the above has been completed, the actual simulation of external conversion proceeds in the user-provided subroutine G U S T E P , called after each tracking step. Appendix D lists a sample G U S T E P modified for modelling external conversion and Compton scattering. The general procedure is as follows: 1. Generate a photon. B y access to the flag I T R T Y P ( G E A N T common G C K I N E ) we can test if we have a photon. It is usually only necessary to convert firstgeneration 7's; second-generation conversions produce electrons of too low energy. Chapter 3. Detector Simulation: Tracking 40 •024 o .021 -Sl2 - .018 ? o - .015 £ 8- c o * 6- - .012 o 3 .009 — liquid h y d r o g e n B o tr h .006 .o c 5 - .003 % 0. .000 2 — I 1 50 150 100 Gamma 200 250 300 energy (MeV) Figure 8: Total cross section for pair production (barn/atom) as calculated using G E A N T parameterized fits of cross section data . 40 60 Gamma energy (MeV) r 80 Figure 9: Total cross section for Compton scattering (barn/atom) as calculated using G E A N T parameterized fits of cross section data . Chapter 3. Detector Simulation: 41 Tracking 2. Decide where the photon is to convert, (random number) < p\. 3. Step the photon through the detector to the chosen detector part. When the 7 reaches this volume, force conversion (by setting the step length to the next pair production to 0) inside the volume, and calculate the probability for the process to have occurred at that point. The relevant quantity is again the interaction length. G E A N T will calculate this energy-dependent quantity using internal lookup tables. The example code found in the appendix uses source directly lifted from G E A N T . In the case of very thin volumes (when it is not necessary to know exactly where the 7 converted) the desired probability is simply the distance to the next boundary, S N E X T (found using G N E X T ) , divided by the interaction length. This gives the number of interaction lengths, which is equal to the probability for one interaction when S N E X T < L. In cases where it is important to model different locations for conversion inside the medium, the step length S T E P must be changed to some random distance between 0 and S N E X T . The weight of the event is then S T E P / L . G E A N T will then step the photon to the location of next interaction and generate the interaction. Unfortunately, G E A N T does not allow for the user to change the step length. S T E P is calculated by G E A N T itself in the subroutine G T G A M A , which is the track-stepping subroutine for photons. One solution is to override the G E A N T subroutine with one in which the user is allowed to change S T E P , see appendix D for a listing. Another solution is to allow G E A N T to set the step length (to S N E X T ) , and hence to force the interaction at the volume boundary, but then to modify the event vertex. G T G A M A will call the subroutines G P A I R G or G C O M P , depending on whether the user forces an external conversion or a Compton scatter. G P A I R G generates an e e~ pair by sampling the Coulomb-corrected Bethe-Heitler cross-section, and G C O M P + Chapter 3. Detector Simulation: Tracking 42 samples the scattered photon energy from the Klein-Nishina formula. See the G E A N T user's guide [14] and the references therein for further details. In general, in order to calculate the interaction lengths, G E A N T makes use of fits to data, the parameters of which have been stored in internal tables. The user's guide gives more details on these fits as well as outlining their ranges of credibility. Chapter 4 Conclusions In this document, we have attempted to present a basic introduction to the techniques involved in writing computer simulations for particle physics experiments concerning rare decays, in which case it is important to understand the functioning of the experimental equipment to a very high degree of accuracy. The reader should now have a feeling for the general procedure and the extent of the work involved, which is not be taken lightly. We give a final brief summary. In order to write a simulation, the programmer must first of all be familiar with the physics of the experiment; enough of the theory should be understood to be able to understand quantitatively the behaviour of the cross section functions, in order to test the sampling routines. The subsequent effects of the kinematics (Lorentz boosts) must be at least qualitatively grasped. Writing the detector simulation usually requires a great deal more effort. Before even beginning, the exact configuration of the detector should be clearly established. The principles of operation and output of each of the detecting elements must be understood, so that the hits and their digitization can be properly handled. The format in which the experimental data will be stored must also be known, in the case that the analysis of the experiment is to be modelled (during the final analysis), so that the same analysis programs can be run on both sets; the same selection criteria and cuts applied to the real and simulated data. We note here that work on a Monte Carlo simulation is never complete. Once 43 Chapter 4. Conclusions 44 analysis of the experimental data begins, the simulation will be tuned and re-tuned as resolutions of various pieces of electronics are matched, as new theoretical results change the event generation, as new software becomes available, etc. During the analysis there is a constant interplay between simulation and real data; the Monte Carlo results are used to understand reality, giving clues as to tuning experimental cuts to eliminate as much of the background as possible. When the data sample is as clean as possible, the roles are reversed and data becomes the final arbiter; if the simulation does not match the data, it is a clue that something is not being properly modelled, and the code must be improved. In appendix E the reader can see some (preliminary) results of a Monte Carlo simulation. A n experiment was done using a beam of 7r 's incident on a liquid hydrogen tar_ get, as described in section 2.3. The real data, consisting essentially of the 4-momenta of any electrons found in each event (from any of the processes 2.9 through 2.13), have been skimmed using rough cuts in order to isolate a sample of TV° —> e + e~ + 7 events. + Geometric cuts on the event vertex help to eliminate much of the external conversion background, while kinematic constraints are used to separate the desired decay from the other reactions. The plots 12 and 13 show simulated and real distributions of x and y (see equations 2.7 and 2.8); the agreement is quite good, indicating that the understanding of both the physics involved and the response of the detector are well understood. Bibliography [1] C. Cohen-Tannoudji, B . D i u , F . Laloe, Quantum Mechanics, J . Wiley and Sons, 1977 [2] F . Halzen, A . D . Martin, Quarks and Leptons: an Introductory Course in Modern Particle Physics, J . Wiley and Sons, 1984 [3] N . M . K r o l l , W . Wada: Phys. Rev. 98, 1355 (1955) [4] T. Miyazaki, E . Takasugi: Phys. Rev. D 8, 2051 (1973) [5] D . W . Joseph: Nuovo Cimento 16, 2021 (1960) [6] F . James, Monte Carlo Phase Space, C E R N publication 68-15 (1968) [7] F . James: Rep. Prog. Phys. 43, 1145 (1980) [8] D I V O N 4 (long writeup), C E R N Computing Center Program Library, D151 [9] J . H . Friedman, M . H . Wright: A C M Trans. Math. Soft. 7, 76 (1981) [10] R. Kleiss, W . J . Stirling, S. D. Ellis: Comp. Phys. Comm. 40, 359 (1986) [11] F O W L (long writeup), C E R N Computing Center Program Library (6000 Series), W505 [12] R. Hagedorn, Relativistic Kinematics, W . A . Benjamin, Inc., 1963, 1973 [13] M . M . Block: Phys. Rev. 101, 796 (1956) [14] G E A N T 3 User's Guide, C E R N Data Handling Division, D D / E E / 8 4 - 1 (1986) 45 Bibliography 46 [15] Particle Data Group, Review of Particle also the associated Particle Properties Properties: Data Booklet Phys. Let. 1 7 0 B (1986) and Appendix A D I V O N coded example The following shows an example of the use of D I V O N . Please note well the use of REAL*8 and REAL*4. The subroutine DIVQNDZ is designed to be called several times; the first time partitioning, integrating, and sampling; the subsequent times sampling only. The variables in the common CUTS are read in in the calling program and can hence be changed without having to recompile and link. Various D I V O N steering switches are also read in. Attention is also directed to the use of ETA in DFUN, which is an example of a variable limit of integration (Y from -ETA to +ETA), and also to EMIN, an example of a further cut on the integration space. 47 Appendix A. DIVON coded example SUBROUTINE DIVONDZ(XVECTOR) Q*********************************************************** IMPLICIT REAL*8 ( A - H , 0 - Z ) REAL*4 SPRDMAX, X V E C T 0 R ( 2 ) , X M I N ( 2 ) , XMAX(2) COMMON / QUADRE / IDEG COMMON / CONSTANTS / P I , RR, GP, GN, ASLOPE COMMON / CUTS / XMINIMUM, XMAXIMUM, EMIN, PETMIN, PPTMIN, ANGMIN, ANGMAX LOGICAL INIT DATA INIT / . F A L S E . / , N / 2 / IF (INIT) C initialize DIVON i n t e g r a t i o n XMIN(1) XMAX(l) C limits GOTO 10 parameters, limits on X = XMINIMUM = XMAXIMUM on Y XMIN(2) = - S Q R T ( 1 . - RR) XMAX(2) - -XMIN(2) C control 101 parameters OPEN (UNIT=10, READONLY, SHARED, READ ( 1 0 . 101) SPRDMAX, MAXPTS F O R M A T ( F 1 0 . 0 , 110) C partition integration space CALL PARTN(N, XMIN, XMAX, C estimate 102 STATUS='OLD') SPRDMAX, MAXPTS) integral READ ( 1 0 , 102) JDEG, NPT FORMAT(2I10) CALL INTGRL(N, JDEG, NPT, ANS, ERROR) CLOSE(UNIT=10) INIT = . T R U E . C f i l l a r r a y XVECTOR ( v a l u e s o f X , Y ) w i t h C v a l u e s a c c o r d i n g t o f u n c t i o n DFUN 10 CALL RANGEN(N, RETURN END XVECTOR) random X , Y Appendix A. DIVON coded example DOUBLE PRECISION FUNCTION DFUN(N, XY) c*********************************************************** C j o i n t p r o b a b i l i t y f n . for X and Y, from K r o l l and Uada IMPLICIT REAL*8 (A-H, 0-2) REAL*8 MP, MN, MPO, MPM, ME, MAXEN REAL*8 XY(2) COMMON / USERMASSES / MP, MN, MPO, MPM, ME COMMON / CUTS / XMINIMUM, XMAXIMUM, EMIN, PETMIN, PPTMIN, ANGMIN, ANGMAX COMMON / CONSTANTS / P I , RR, GP, GN, ASLOPE DATA CONST / 5.80857D2 / X = XY(1) Y = XY(2) IF IF ETA IF (X . L T . RR) GOTO 100 ! redundant but safe (X . L T . XMIN .OR. X . G T . XMAX) GOTO 100 = SQRTC1. - RR/X) (ABS(Y) . G T . ETA) GOTO 100 EPAIR PPAIR EP0 = EE0 = = MPO/2. * (1. + X) = MPO/2. * (1. - X) (EPAIR + PPAIR*Y) / 2. EPAIR - EP0 IF (EP0 . L T . EMIN .OR. EE0 . L T . EMIN) GOTO 100 100 FORM0 = 1.D0 + 2.D0*ASLOPE*X DFUN = CONST * (1.D0 - X)**3 * (1.D0 + RR/X + Y*Y) * (1.D0/X) * FORM0 RETURN DFUN = 0.0D0 RETURN END Appendix B G E A N T detector d e f i n i t i o n The following code is an example of the definition of a particular detector, shown in figure 10. This subroutine is called by U G I N I T . The materials and mixtures to be used are defined first. Then these materials are used to define the tracking media. In defining materials, the user supplies a material number, the first number in the GSMATE call. This number is used in setting the tracking media, and appear as the third number is the GSTMED calls. The names supplied are not used in G E A N T and are for user convienience. For example, in S I N D R U M , the material S C I N T is used to define the tracking media P O L Y S T Y R E N E and P 0 L Y H 0 D 0 . P 0 L Y H 0 D 0 is used in the definition of the hodoscopes. The geometry parameters (number of wires per wire chamber, relative position of the chambers, chamber wall thickness, etc.) are read in from a user-defined data bank EGEO, as some of these are subject to change as the detector is taken apart, repaired, and put back together over the course of the experimental run. Note that all dimensions are in cm. These parameters are used to define and position the volumes of S I N D R U M . Notice that this simulation does not make use of the G E A N T hit and digitization facilities, using user-written routines instead; there are no calls to G S D E T . 50 Appendix B. GEANT detector definition Figure 10: The S I N D R U M detector. t 51 Appendix B. GEANT detector definition SUBROUTINE GEOM C Define SINDRUM geometry set up and materials £*********************************************************** DIMENSION AMYLAR(3),ZMYLAR(3),WMYLAR(3) DIMENSION AMAKRO(3),ZMAKROC3),WMAKROC3) DIMENSION ASC1NT(2),ZSCINT(2),USCINT(2) DIMENSION AKAPA(4),ZKAPA(4),UKAPA(4) DIMENSION PAR(8) CHARACTERS CHANAM(7),GASNAM(7),UIRNAM(7),ACTNAM(7) INCLUDE 'SIND$GEANTCOMM(GCONST)' INCLUDE 'SINDSGEANTCOMM(GCFLAG)' INCLUDE 'SINDSGEANTCOMM(GCBANK)' INCLUDE 'INCLUDESDIRrBCS.FOR' INCLUDE 'INCLUDE$DIR:CHODO.FOR' INCLUDE 'INCLUDE$DIR:CCHAM.FOR' INCLUDE 'INCLUDE$DIR:CCATH.FOR' COMMON/FIELD/FIELDM COMMON/CEXPNU/NUMEXP,NSTEPS,NACTUAL,NRUN C data for s c i n t i l l a t o r and "kapacell" compound DATA DATA DATA DATA DATA DATA DATA DATA DATA DATA ASCI NT/12.011,1.008/AKAPA/12.,1.,14.,16./ ZSCINT/6.,1./ZKAPA/6..1..7.,8./ USCINT/1.,1.1/WKAPA/19.,14.,2.,4./ AMYLAR/12..1..16./AMAKRO/12.,1.,16./ ZMYLAR/6.,1.,8./ZMAKR0/6.,1.,8./ WMYLAR/5.,4.,2./WMAKRO/3..8.,2./ CHANAM/'CHA1','CHA2','CHA3','CHA4','CHA5','CHA6','CHA7'/ GASNAM/'ARG1','ARG2','ARG3','ARG4','ARG5','ARG6','ARG7'/ ACTNAM/'ACT1','ACT2','ACT3','ACT4','ACT5','ACT6','ACT7'/ UIRNAM/'WIR1','wIR2','UIR3','WIR4','WIR5','wIR6','wIR7'/ C read i n SINDRUM geometry from EGEO bank CALL READ_EGEO C Define materials CALL CALL CALL CALL CALL CALL 1.01, 1.,0.071 . 8 6 5 . . 7 9 0 . , 0 , 0 ) 4.00, 2.,0.000178,1075120.,680.E3,0,0) 26.98,13.,2.7 , 8.9,37.2,0,0) 55.85,26.,7.87 ,1.76,17.1,0,0) 63.54,29.,8.96 .1.43.14.8,0,0) 14.61,7.3,0.001205,30423.,67500.,0,0) GSMATE( 1,'HYDROGENS GSMATE( 3,'HELIUMS GSMATEC 9,'ALUMINIUMS GSMATE(10,'I RONS GSMATEC11,'COPPERS GSMATE(15,'AIR$ C f i x up GEANT bug with vacuum CALL GSMATEC16,'VACUUMS CALL CALL CALL CALL CALL ', 4.00, 2.,0.00000178,1.E8,1.E8,0,0) GSMATE(17,'ARGONDENSS', 39.948,18.,0.00178,16500.,0.,0,0) GSMIXTC20,'SCINTS ',ASCI NT,ZSCINT,1.032,-2,WSCI NT) GSMIXT(21,'KAPACELLS ',AKAPA,ZKAPA,0.15,-4.WKAPA) GSMIXT(22,'MYLARS ',AMYLAR,ZMYLAR,1.39,-3,WMYLAR) GSMIXT(23,'MAKROLONS ',AMAKRO,ZMAKRO,1.18,-3,WMAKRO) C Defines USER tracking media parameters FIELDM I FIELD TMAXFD STMIN = 3.313 = 3 = 10.0 = 0.0000001 DMAXMS = DEEMAX = EPSIL = * 0.5 0.2 0.01 !FIELD STRENGTH !FIELD CONSTANT !MAX. ANGLE DUE !MIN. STEP DUE IN kG ALONG AXIS 3 TO FIELD (DEGREES) TO MULT. SCATTERING (CM) !MAX. DISPLACEMENT DUE TO MULT. SCATTERING (CM) !MAX. FRACTIONAL ENERGY LOSS !TRACKING PRECISION (CM) CALL GSTMED( 1,'DEFAULT MEDIUM HELIUMS',3,0,1 FIELD,FIELDM, TMAXFD,DMAXMS,DEEMAX,EPSIL,STMIN, 0,0) DMAXMS = DEEMAX = EPSIL = 0.01 0.003 0.0003 !MAX. DISPLACEMENT DUE TO MULT. SCATTERING (CM) !MAX. FRACTIONAL ENERGY LOSS ITRACKING PRECISION (CM) Appendix B. * * * * * * * * * * definition CALL GSTMEDC 2,'ARGONS' ,17,0,1 FIELD,FIELDM, TMAX FD,DMAXMS,DE EMAX,EPSIL,STMIN,0,0) CALL GSTMEDC 3,'POLYSTYRENES' ,20,0,1 FIELD,FIELDM, TMAXFD,DMAXMS,DEEMAX.EPSIL.STMIN,0,0) CALL GSTMEDC 7,'KAPACELLS' ,21,0,1 FIELD,FIELDM, TMAXFD,DMAXMS,DEEMAX.EPSIL,STMIN,0,0) CALL GSTMEDC 8,'HYDROGENS' , 1,0,1 FIELD,FIELDM, TMAXFD,DMAXMS,DEEMAX,EPS IL,STMIN,0,0) CALL GSTMEDC22,'MYLARS ' ,22,0,1 FIELD,FIELDM, TMAXFD.DMAXMS,DEEMAX,EPS IL,STMIN,0,0) CALL GSTMED(23,'MAKROLONS' ,23,0,1 FIELD,FIELDM, TMAXFD,DMAXMS,DEEMAX,EPS IL.STMIN,0,0) DMAXMS = DEEMAX = EPSIL = * GEANT detector 0.5 0.2 0.01 !MAX. DISPLACEMENT DUE TO MULT. SCATTERING CCM) !MAX. FRACTIONAL ENERGY LOSS ITRACK1NG PRECISION CCM) CALL GSTMEDC16,'VACUUMS ' ,16,0,1 FIELD,FIELDM, TMAXFD,DMAXMS,DEEMAX,EPSIL,STMIN,0,0) CALL GSTMEDC10,'IRONS ' ,10,0,1 FIELD,FIELDM, TMAXFD,DMAXMS,DEEMAX,EPSIL,STMIN,0,0) CALL GSTMEDC A,'ALUMINIUMS' , 9,0,1 FIELD,FIELDM, TMAXFD.DMAXMS,DEEMAX,EPS IL.STMIN,0,0) CALL GSTMEDC 5,'COPPERS' ,11,0,1FI ELD,FIELDM, TMAXFD,DMAXMS,DEEMAX,EPS IL.STMIN,0,0) CALL GSTMEDC 6,'POLYHODOS' ,20,0,IFIELD,FIELDM, TMAXFD,DMAXMS,DEEMAX,EPS IL.STMIN,0,0) C Define volumes using numbers from EGEO C mother volume = master reference system CMARS) PARCD=0. PARC2)=41. PARC3)=70.0 CALL GSV0LU('SIND','TUBE',1,PAR,3,IV0L ) C chambers DO I=1,NCHA PARCl)=RINCI)-DWALLCI) PARC2)=RINCI)+DWALLCI)+GAPCI) PARC3)=CZUPCI)-ZDOUNCI))*0.5 CALL GSVOLUCCHANAMCI),'TUBE',7,PAR,3,1 VOL) PARCD=RINCI) PARC2)=RINCI)+GAPCI) CALL GSVOLUCGASNAMCI),'TUBE',2,PAR,3,IVOL) PARC1)=RINCI)+(GAPCI)-ACTIVGAPCI))*0.5 PARC2)=RINCI)+(GAPCI)+ACTIVGAPCI))*0.5 CALL GSVOLUCACTNAMCI),'TUBE',2,PAR,3,1 VOL) ENDDO C hodoscope: PARC1)=HODRAD PARC2)=HODRAD+HODTHI PAR(3)=H0DLEN/2. CALL GSVOLUC'HODO','TUBE',6,PAR,3,1 VOL) C aluminium c y l i n d e r : PARC1)=35.6 PARC2)=38.4 PAR(3)=66.0 CALL GSV0LUC'ALUM','TUBE',A,PAR,3,IV0L ) C magnet c o i I : PARC1)=38.4 PARC2)=40.7 PARC3)=66.0 CALL GSV0LUC'C0IL','TUBE',5,PAR,3,IV0L ) 53 Appendix B. GEANT detector definition C P o s i t i o n volumes within SIND (MASTER REFERENCE) CALL TARDEF IDEFINES AND POSITIONS TARGET DO I=1,NCHA 20=(ZUP(I)+ZDOUN(I))/2. CALL GSPOS(CHANAM(!),1,'SIND',0.,0.,Z0,0,'ONLY') CALL GSPOS(GASNAM(I),1,CHANAM(I),0.,0.,0.,0,'ONLY') CALL GSPOS(ACTNAM(I),1,GASNAM(I),0.,0.,0.,0,'ONLY') ENDDO CALL G S P O S ( ' H O D O ' , 1 , ' S I N D ' , 0 . , 0 . , 0 . , 0 , ' O N L Y ' ) CALL G S P O S ( ' A L U M ' , 1 , ' S I N D ' , 0 . , 0 . , 0 . , 0 , ' O N L Y ' ) CALL G S P O S ( ' C O I L ' , 1 , ' S I N D ' , 0 . , 0 . , 0 . , 0 , ' O N L Y ' ) C Divide chambers and hodoscope DO I=1,NCHA C0=-I80./NWIRE(I)+PHIO(I) CALL GSDVN2(WIRNAM(I),ACTNAM(I),NWIRE(I),2,C0,2) ENDDO C0=-180./NSECT+HODPHI0 CALL GSDVN2('SECT','HODO',NSECT,2,C0,6) C order 1st-generation daughters of SIND for optimal tracking CALL GSORD('SIND',4) CALL GSORD('VACU',4) C Close geometry banks, ( o b l i g a t o r y system routine) CALL GGCLOS RETURN END SUBROUTINE READ EGEO £*********************************************************** C Get SINDRUM geometry from EGEO bank £*********************************************************** INCLUDE 'INCLUDE$DIR:BCS.FOR' INCLUDE 'INCLUDESDIR:CHODO.FOR' INCLUDE 'INCLUDESDIR:CCHAM.FOR' INCLUDE 'INCLUDE$DIR:CCATH.FOR' C chamber r o t a t i o n s for various runs 1 2 3 4 5 6 7 8 REAL DPHINEU(5, 8) INTEGER NEUGEO(8) DATA NEWGEO / 0, 92, 160, 171, 208, 225, 265, 270 / DATA DPHINEW / 0.41425, 1.69182, 0.1199. 0.28887. 0.17578. 0.45825, 1.72272, -0.01852, 0.22555, 0.17578, -0.259 , 1.71401, 0.04749, 0.09988, 0.17578, 2.52689, 1.71401, 0.04749, 0.09988, 0.17578, 1.404 , 1.68237, 0.02747, 0.0992, 0.17578. 1.9934 , 1.71338, 0.04483, 0.09966, 0.17578, 2.21378, 1.97664, 0.20499, 0.39048, 0.17578, 2.26254, 2.00748, 0.46625, 0.53028, 0.17578 / C f i n d out geometry for t h i s run 10 1000 DO 111=1,8 IF(ABS(NRUN).LE.NEWGEO(IU)) GOTO 10 END DO IU=IU-1 WRITE(6, 1000) IU FORMATC ROTATION SET USED:', 1 2 , / , / ) Appendix B. GEANT C C C C C C detector definition f i l l geometry commons with data from bank EGEO (for chambers and hodoscope) ATTENTION: only one hodoscope supported; angle between anode and inner cathode is assumed to be zero as well as z - p o s i t i o n of middle of chamber for chambers with cathode readout. IP=IDATA(NAMIND('EGEO')) L0=IDATA(IP+1) JH0D=IDATA(IP+2) LTRHOD=IDATA(IP+3) NCHA=IDATA(IP+4) LTRCH=IDATA(IP+5) IF(JH0D.GT.1)ST0P 'ONLY ONE HODOSCOPE ALLOWED' IF(NCHA.GT.7)ST0P 'ONLY 7 CHAMBERS ALLOWED' IP=IP+L0 H0DRAD=RDATA( IP+D/10 NSECT=IDATA(IP+2) HODPHI0=RDATA(IP+3)*RADDEG HODLEN=RDATA(IP+4)/10 HODTHI=RDATA(IP+5)/10 IP=IP+LTRHOD ! HODOSCOPE GEOMETRY DO ICH=1,NCHA !CHAMBER GEOMETRY RIN(ICH)=(RDATA(IP+1)-RDATA(IP+5))/10 GAP(ICH)=2*RDATA(IP+5)/10 ACTIVGAP(ICH)=0.15 DWALL(ICH)=0.31*RDATA(IP+6) ZDOWN(ICH)=(RDATA(IP+14)-RDATA(IP+4)/2)/10 ZUP(ICH)=(RDATA(IP+14)+RDATA(IP+4)/2)/10 NWIRE(ICH)=IDATA(IP+2) PHIO(ICH)=DPHINEW(ICH, IU) BOTHSIDES(ICH)=.FALSE. IF(IDATA(IP+7).EQ.'N0')THEN !INNER CATHODES NSTRIPS(ICH,1)=0 ELSE IHIGFIR=IDATA(IP+9)/65536 ILOWFIR=IAND(IDATA(IP+9),65535) ILOWLAS=IAND(IDATA(IP+10),65535) IF(IHIGFIR.NE.0)BOTHSIDES(ICH)=.TRUE. NSTRIPS(ICH,1)=IL0WLAS-IL0UFIR+1 END IF IF(IDATA(IP+8).EQ.'N0')THEN NSTRIPS(ICH,2)=0 ELSE IHIGFIR=IDATA(IP+11)/65536 ILOWFIR=IAND(IDATA(IP+11),65535) ILOWLAS=IAND(IDATA(IP+12),65535) IF(IHIGFIR.NE.0)BOTHSIDES(ICH)=.TRUE. NSTRIPS(ICH,2)=ILOWLAS-ILOWFIR+1 ENDIF IP=IP+LTRCH ENDDO RETURN END !OUTER CATHODES Appendix C G E A N T t r a c k i n g code In the case of the S I N D R U M simulation, we do not make use of G E A N T ' s digitization routines. Hits are recorded in a different bank structure, whose management routines begin with ' B ' . These banks are later used in the digitization, and are used for the online data also. In GUTREV, counters for the number of track steps and the number of hits per track are initialized, as well as the logical variable HODOREACHED. The hit bank-structure is cleared. Subroutine GUSTEP aborts events with too many tracking steps, stores only those secondary particles which have a chance of reaching sensitive detector elements (those not generated in the hodoscopes or magnet coil surrounding S I N D R U M ) , and stores hits in various arrays. The wire number hit is stored, also the energy deposited in the contact. The G E A N T array NAMES (here equivalenced to the character array TEXT) gives the name of the last tested volume. Unless a wire has been hit, no hit information is stored. For a track just entering a volume ( G E A N T variable INWVOL = 1), the chamber and wire numbers are stored, as well as the z-component of the momentum upon entering the chamber. If the track is entering a hodoscope sector, the sector number is stored. When the track leaves the volume, the z-momentum is again stored, and the energy deposited in the volume is accumulated (DESTEP is the total energy lost in the current step). 56 Appendix C. GEANT tracking code SUBROUTINE GUTREV c*********************************************************** C User routine to control tracking of one event C C a l l e d by GRUN C*********************************************************** INCLUDE 'SIND$GEANTCOMM(GCNUM)' INCLUDE 'SIND$GEANTCOMM(GCFLAG)' INCLUDE 'SIND$GEANTCOMM(GCKINE)' INCLUDE 'INCLUDE$DIR:BCS.FOR' INCLUDE 'INCLUDE$DIR:CHITS.FOR' INCLUDE 'INCLUDE$DIR:CCHAM.FOR' LOGICAL HODOREACHED COMMON/CEXPNU/NUMEXP,NSTEPS.NACTUAL,NRUN C0MM0N/CSTEER/H0D0REACHED(10) NACTUAL=0 NUMTRA=NTRACK DO I=1,NUHTRA NHITS(I)=0 HODOREACHED(I)=.FALSE. ENDDO !number of tracking s CALL BDROP(IDATA,'E') CALL BGARBCIDATA) !clear out h i t bank CALL GTREVE !do the tracking RETURN END Appendix C. GEANT tracking code SUBROUTINE GUSTEP C*********************************************************** C User routine c a l l e d at the end of each tracking step C INUVOL i s d i f f e r e n t from 0 when the track has reached C a volume boundary C ISTOP i s d i f f e r e n t from 0 i f the track has stopped C*********************************************************** LOGICAL INCHAMBER,INHODOSKO CHARACTERS TEXT(15) EQUIVALENCE (NAMESC1),TEXT(1)) INCLUDE 'SIND$GEANTCOMM(GCKING)' INCLUDE 'SIND$GEANTCOMM(GCKINE)' INCLUDE 'SIND$GEANTCOMM(GCTRAK)' INCLUDE 'SIND$GEANTCOMM(GCFLAG)' INCLUDE SIND$GEANTCOMM(GCTMED)' INCLUDE 'SIND$GEANTCOMM(GCNUM)' INCLUDE 'SINDSGEANTCOMM(GCVOLU)' INCLUDE 'SINDSGEANTCOMM(GCBANK)' INCLUDE 'INCLUDE$DIR:CCHAM.FOR' INCLUDE 'INCLUDE$DIR:CHITS.FOR' COMMON/CEXPNU/NUHEXP,NSTEPS,NACTUAL,NRUN LOGICAL HODOREACHED COMMON/CSTEER/HODOREACHED(10) ; C count number of tracking steps - abort event i f too many C steps 900 NACTUAL=NACTUAL+1 IFCNACTUAL .GT. NSTEPS)THEN WRITE(6,900) IDEVT FORMAT(2X, 15, ' EVENT ABORTED - TOO MANY TRACKING STEPS') IEOTRI=1 RETURN END IF C store secondaries into stack - - only those o r i g i n a t i n g C inside target or chambers IF(NGKINE.GT.O) THEN IF(TEXTCNLEVEL) . N E . 'COIL' .AND. TEXT(NLEVEL) . N E . 'ALUM' .AND. TEXT(NLEVEL) .NE. 'SECT') CALL GSKING(O) ENDIF C produce h i t data INCHAMBER=TEXT(NLEVEL)(:3).EQ.'WIR' INHODOSKO=TEXT(NLEVEL).EQ.'SECT' IF(INHODOSKO)HODOREACHED(ITRA)=.TRUE. C are we in a s e n s i t i v e volume—chambers or hodos only? IF(.NOT.(INCHAMBER.OR.INHODOSKO))RETURN C 100 h i t s only IF(NHITS(ITRA).GE.100)RETURN C only interested i n e+ or e- p a r t i c l e s in gas IF(INCHAMBER.AND.ITRTYP.NE.2)RETURN C just entering volume? IF(INUVOL.EQ.1)THEN C store number of h i t s in t h i s track, Pz at entry point C (P from GEANT bank, wire # h i t , and chamber/hodo #. NHITS(ITRA)=NHITS(ITRA)+1 NH=NHITS(ITRA) ZENTRY(ITRA,NH)=VECT(3) NUMWIR(ITRA,NH)=NUMBER(NLEVEL)-1 IWIRES START AT 0 EDEPO(ITRA,NH)=0. IF(INHODOSKO)THEN NUMCHACITRA,NH)=NCHA+1 ELSE NUMCHA(ITRA,NH)=ICHAR(TEXT(NLEVEL)(4:4))-ICHAR('0') ENDIF ELSE Appendix C. GEANT tracking code C or e x i t i n g volume? IF(INUVOL.GE.2.OR.ISTOP.NE.0)THEN C then store Pz at e x i t p o i n t , add up E dep. inside volume C DESTEP from GEANT bank NH=NHITS(ITRA) ZEXIT(1TRA,NH)=VECT(3) EDEPOCITRA,NH)=EDEPO(ITRA,NH)+DESTEP C or are we s t i l l inside the volume? --then just add up E ELSE EDEPOCITRA,NH)=EDEPO(ITRA,NH)+DESTEP ENDIF END IF RETURN END Appendix D E x t e r n a l C o n v e r s i o n code Here we list sample code for external conversion and Compton scattering. Only primary gammas are converted/scattered, these may be from either of the two processes (D.16) 7T° -> 7 + 7 7T° -> e + e-+ + 7 (D.17) Subroutine GUTRAK is called at the start of each track. NG is the limit for the number of gammas to track; if NG > 2 we do not have a primary 7 anymore. For D.17 there is only one primary photon, so we stop after converting it. Subroutine GUSTEP is called after each tracking step. This is where conversion/scattering occurs. Array NYTPE is read in previously, and specifies which process to model for each of the NG gammas; NTYPE = 1 for Compton scattering, NTYPE = 2 for conversion. Array NWHERE gives the location for converting/scattering each gamma. NEWEVENT and DALITZ are set every time a new event is generated in GUKINE. A substitute version of the G E A N T subroutine GTGAMA, the photon track-stepping routine, is also given. It is modified to allow the user to redefine the track step, by setting the flag EXTCON and thereby forcing G E A N T to override STEP with the userprovided USER_STEP. 60 Appendix D. External Conversion code Q*********************************************************** SUBROUTINE GUTRAK C User routine to control tracking at track INCLUDE 'lNCLUDE$dir:CSTEER.FOR' INCLUDE 'sind$geantcomm(gcflag)' INCLUDE 'sind$geantcomm(gckine)' INCLUDE 'INCLUDESdi r:CKINE2.FOR' C contains NTYPE, NG, NWHERE INCLUDE 'DB1:[CEW.SDM]CTYPE.FOR' C i s t h i s a photon track? IF(ITRTYP.EQ.1)THEN IF(NEWEVENT)THEN NG = 1 I F(DALITZ)NG = 2 WEIGHT = 0. WT(1) = 0. WT(2) = 0. ELSE NG = NG+1 END IF NEWEVENT = .FALSE. END IF CALL GTRACK RETURN END level. (A13dVS'lX3NS'l33A)lX3N3 11V3 1X3NS Ajepunoq j x s u 01 s o u e a s i p 3 I 0 1 03((3N)3a3HnrT3N"3a3HrU)dI 1 = 3a3HnI(;lVH3/"03 O3A31N)lX3i)dI £ = 3a3Hni(,NinA, D3"(T3A31N)XX31)dl Z = 3d3Hni(i1f1i,-03'(£:)(13A31N)iX31)dI I = 3a3HnI(,HVl, 03-(£:)O3A31N)lX31)dI 0 = SaSHttl - - - 11 136J0J '(3H3HMN) 3 U01SJ3AUCO 33J0J. 0} JUBM 3M 3J3UM JOU 3J,3M MOU 3JB 0 a n a j a i j M s i l a j 3a3HrU ' u o i s j a A u o o a o j o j oi a J a u M a p p a p 3 •3nbl' = N031X3 1NV33 JOJ- 6 e i j . uo.isjaAUOo i s u j a j x a } a s o s 'saA 3 I 0 1 03(0'03'(ON)3dAiN)dI eimis6 Ajeuiud e ION i 1 o i 03<Z'13'3N)dI I 01 03(r3N dAiHira0T3N 10ArlNI)dI - - iauinioA M3u e 6u_u3iU3 j s n f 11 s i pus BIUUJBS B SABU SM oa o uojduio3 j o uoijonpojd J i e a aojoj 3 d l 0N3 I = I8103I / i s o a v i s a i s i n o a s i a v i s iN3A3/'*3dAi N3Hl(,aVl,-3N'(£:)(13A31N)lX31 aNW •0"O3 3N31S aNV£ O3'dAiaiI)dI - - - - 1S6JEJ ei^i u i p s } j B ) s u s a q s e q JUSAG am >(33qo 3 dI0N3 Nanisa 1 = iaiosi (,Sd31S ONIXOVMl ANVW 001 - 0 3 i a o a V < 'X£'9I'/1N3A3 /HVWHOd N3Hi(Sd31SN "13" "IVni3VN)dI l+1Vni3VN = 1VP13VN 006 s d 3 J S AUBUI 001 3 - s d s j s 6 u i 5 | 3 e j j j.o jsquinu j u n c o 3 j.i JU3A3 u o q s d31S~il3Sn 'N031X3 / N031X3 / N0WW03 (Ol)a3H3V3»OaOH/a331S3/NOWW03 N031X3 'CBHDtfSaoaOH 1V3I301 1Vni3VN'Sd31SN'dX3WnN/nNdX33/N0WW03 ,aOd 3dA13[WaS M33] '-IBQ, 3am3NI iHOd"Z3NI>3:aia$30m3NI / 3 a m 3 N I /aod siiH3:aia$3ani3Ni; aamoNi /aod wvH33:aia$3arn3Ni, gamoNi ,(rnoA3D)wwo3iNV33$aNiS/ samoNi /(WnN33)WWO31NV30$aNIS; samoNi (a3W133)WW031NV3D$aNISi s a r m m i(3V1d33)WW031NV3D$aNIS( 3am3NI ;(SAHd03)WW031NV33$aNISi SuTTONI /(>Vai33)WW031NV3D$aNISi 3afll3NI i(3NIX33)WW031NV33$aNIS/ 3 a m 3 N I /(3NI)f33)WW031NV33$aNIS/ SamONI j(30ir3D)WW031NV33$aNIS/ 3 a m 3 N I /(>INVa33)WW031NV33$aNIS/ samoNi - - - - ( ((L)1X31'(L)S3WVN) 33N31VAI003 (51)1X31 t7»M313VaVH3 oxsoaoHNiasawvHONi "iV3iooi ***********************************************************3 * p s d d o j s s e q yoeji j.t 0 WOJJ. j u a j a j j i p s j dOlSI 3 » Ajepunoq adinioA e 3 » P3U0B3J SBU, >|3BJ1 3U,} U3L|M 0 U10J}. 1U3J3J.JIP SI 10AMNI » d s j s 6u_i>|3BJi ipsa j.o pua ai|i IB p a n s a a u u i i o j j s s n 3 ***********************************************************3 d3isn3 B N i i n o a s n s ***********************************************************3 29 9poo UOISI9AUOQ jvujQixg (j xxpusddy 3 Appendix D. External Conversion code C Update pointers temporarily for new medium (geant hasn't C done i t yet) JPAIR-XPAIR, JCOMP-XCOMP KCOMP = LQCJMA-8) KPAIR = LQ(JMA-10) C i n t e r a c t i o n length, code p l a g i a r i z e d from GTGAMA GEKRT1 = 1.-GEKRAT IF(NTYPE(NG).EQ.2)SIGPA = GEKRT1*Q(KPAIR+IEKBIN)+GEKRAT*Q(KPAIR+IEKBIN+1) IF(NTYPE(NG).EQ.1)SIGC0 = GEKRT1*Q(KC0MP+IEKBIN)+GEKRAT*Q(KC0MP+IEKBIN+1) C c a l c u l a t e step s i z e and i n t e r a c t i o n p r o b a b i l i t y for C target: IF(I WHERE. EQ.DTHEN USER STEP = RAN(II)*SNEXT IF(NTYPE(NG).EQ.2)WT(NG) = USER_STEP/SIGPA IF(NTYPE(NG).EQ.1)WT(NG) = USER_STEP/SIGCO C . . . f o r anywhere e l s e : ELSE IF(NTYPE(NG).EQ.2)WT(NG) = SNEXT/SIGPA IF(NTYPE(NG).EQ.1)WT(NG) = SNEXT/SIGCO END IF C Product of weights of two gammas IF(NG.EQ.2)WEIGHT = WT(1)*WT(2) IF(DALIT2)WEIGHT = WT(2) C force conversion. Setting s i . . . and z i n t . . . both forces C distance to next i n t e r a c t i o n to be 1.0E-7 * i n t e r a c t i o n C length. IF(NTYPE(NG).EQ.2)THEN NTYPE(NG) = 0 SLPAIR = SLENG ZINTPA = 1.E-7 END IF IF(NTYPE(NG).EQ.DTHEN NTYPE(NG) = 0 SLCOMP = SLENG ZINTCO = 1.E-7 END IF C store secondaries into stack - - only those o r i g i n a t i n g C inside target or chambers 1 IF(NGKINE.GT.O) THEN IF(TEXTCNLEVEL) .NE. 'COIL' .AND. TEXT(NLEVEL) .NE. 'ALUM' .AND. TEXT(NLEVEL) . N E . 'SECT') CALL GSKING(O) ENDIF C debug requests IF(ISWIT(3).EQ.DTHEN CALL GSXYZ CALL GPCXYZ END IF C produce h i t data INCHAMBER = TEXT(NLEVEL)(:3).EQ.'WIR' INHODOSKO = TEXT(NLEVEL).EQ.'SECT' IF(INHODOSKO.AND.ITRTYP.EQ.2)H0D0REACHED(ITRA) .TRUE. IF(.NOT.(INCHAMBER.OR.INHODOSKO))RETURN {SENSITIVE VOLUME (CHAMBER OR ! HODOSCOPE)? 1ALLOU ONLY 100 HITS IF(NHITS(ITRA).GE.100)RETURN !IGNORE PART. EXCEPT e+- IN GAS IF(INCHAMBER.AND.ITRTYP.NE.2)RETURN !ENTERING VOLUME IF(INWVOL.EQ.DTHEN NHITS(ITRA) = NHITS(ITRA)+1 NH = NHITS(ITRA) ZENTRY(ITRA,NH) = VECT(3) NUMWIR(ITRA,NH) = NUMBER(NLEVEL)-1 IWIRES START AT 0 EDEPO(ITRA,NH) = 0. Appendix D. External Conversion code I F(INHODOSKO)THEN NUMCHA(ITRA.NH) = NCHA+1 ELSE NUMCHA(ITRA,NH) = ICHAR(TEXT(NLEVEL)(4:4))-ICHAR('0') END IF ELSE IF(INWVOL.GE.2.OR.ISTOP.NE.0)THEN "EXITING VOLUME NH = NHITS(ITRA) ZEXIT(ITRA,NH) = VECT(3) EDEPO(ITRA.NH) = EDEPO(ITRA,NH)+DESTEP ELSE !INSIDE VOLUME EDEPOCITRA.NH) = EDEPOCITRA,NH)+DESTEP END I F END I F RETURN END Appendix D. External Conversion code £*********************************************************** C. C c! C. C. C. C. C. C. SUBROUTINE GTGAMA * ****************************************************** * A photon track is extrapolated by one step * * ==>Called by : GTVOL * Authors R.Brun, M.Maire ********* * * Modified by Reena Meijer-Drees 24-06-88 ****************************************************** * * * * * * * PARAMETER (KWBANK=90D0,KWWORK=5200) COMMON/GCBANK/NZEBRA,GVERSN,ZVERSN,IXSTOR,IXDIV,IXCONS,FENDQ(16) + ,LMAIN,LR1,WS(KWBANK) DIMENSION IQ(1).Q(1),LQ(8000) EQUIVALENCE (Q(1),IQ(1),LQ(9)),(LQ(1),LMAIN) COMMON/GCLINK/JDIGI ,JDRAW ,JHEAD ,JHITS ,JKINE ,JMATE ,JPART + ,JROTM ,JRUNG ,JSET ,JSTAK ,JGSTAT,JTMED ,JTRACK,JVERTX + ,JVOLUM,JXYZ COMMON/GCPHYS/IPAIR,SPAIR,SLPAIR,ZINTPA,STEPPA + ,ICOMP,SCOMP,SLCOMP,ZINTCO,STEPCO + ,IPHOT,SPHOT,SLPHOT,ZINTPH,STEPPH + ,IPFIS,SPFIS,SLPFIS,ZINTPF,STEPPF + ,IDRAY,SDRAY,SLDRAY,ZINTDR,STEPDR + , I ANN I,SANNI,SLANNI,ZINTAN,STEPAN + ,IBREM,SBREM,SLBREM,ZINTBR,STEPBR + ,IHADR,SHADR,SLHADR,ZINTHA,STEPHA + ,IMUNU,SMUNU,SLMUNU,ZINTMU,STEPMU + ,IDCAY,SDCAY,SLIFE ,SUMLIF,DPHYS1 + ,ILOSS,SLOSS,SOLOSS,STLOSS,DPHYS2 + ,IMULS,SMULS,SOMULS,STMULS,DPHYS3 + + + COMMON/GCCUTS/CUTGAM,CUTELE,CUTNEU,CUTHAD,CUTMUO,BCUTE,BCUTM ,DCUTE ,DCUTM ,PPCUTM,TOFMAX,GCUTS(5) COMMON/GC FLAG/1 DEBUG, IDEMIN, IDEMAX, ITEST. IDRUN, IDEVT, IEORUN ,IEOTRI,I EVENT,ISWTT(10),IFINIT(20),NEVENT,NRNDM(2) COMMON/GCJLOC/NJLOC(2),JTM,JMA,JLOSS,JPROB,JMIXT,JPHOT,JANNI ,JCOMP,JBREM,JPAIR,JDRAY,JPFIS,JMUNU,JHADR INTEGER + + NJLOC ,JTM,JMA,JLOSS,JPROB,JMIXT,JPHOT,JANNI ,JCOMP,JBREM,JPAIR,JDRAY,JPFIS,JMUNU,JHADR COMMON/GCKINE/HCINE,PKINE(10),ITRA,ISTAK,IVERT,IPART,ITRTYP ,NAPART(5),AMASS,CHARGE,TLIFE,VERTC3),PVERT(4),IPAOLD COMMON/GCKING/KCASE,NGKINE,GKIN(5,100),TOFD(100) INTEGER REAL KCASE,NGKINE GKIN.TOFD COMMON/GCONST/PI.TWOPI ,PIBY2,DEGRAD,RADDEG,CLIGHT ,BIG,EMASS COMMON/GCTMED/NUMED,NATMED(5),ISVOL,IFI ELD,FIELDM,TMAXFD,DMAXMS + ,DEEMAX,EPSIL,STMIN,CFIELD,CMULS,IUPD,ISTPAR,NUMOLD COMMON/GCTRAK/VECT(7),GETOT,GEKIN,VOUT(7),NMEC,LMEC(30),NAMEC(30) + ,NSTEP ,MAXNST,DESTEP,SAFETY,SLENG ,STEP ,SNEXT ,SFI ELD,SNPHYS + ,TOFG ,GEKRAT,IGNEXT,INWVOL,ISTOP ,IDECAD,IEKBIN DATA EPCUT/1.022E-3/ COMMON / EXTCON / EXTCON, USER STEP LOGICAL EXTCON C. C C Is the p a r t i c l e below threshold energy ? IF(GEKIN.LE.CUTGAM)THEN ISTOP=2 DESTEP=DESTEP+GEKIN Appendix D. External Conversion code C C C STEP=EPSIL SLENG=SLENG+EPSIL NMEC=1 LMEC(1)=30 RETURN ENDIF Update local pointers i f medium change IF(IUPD.EQ.O)THEN IUPD=1 JPH0T=LQ(JMA-6) JC0MP=LQ(JMA-8) JPAIR=LQ(JMA-10) JPFIS=LQ(JMA-12) SNPHYS=0. ENDIF C MEC =1 STEP = SNEXT IF(SNPHYS.GT.STEP)GO TO 10 SNPHYS=BIG GEKRT1=1.-GEKRAT C C=======> Get next p a i r production point C IF (GETOT.GT.EPCUT)THEN ZINTPA = ZINTPA - (SLENG-SLPAIR) / STEPPA SLPAIR = SLENG STEPPA=GEKRT1*Q(JPAIR+IEKBIN)+GEKRAT*QCJPAIR+IEKBIN+1) SPAIR = STEPPA * ZINTPA IF(SPAIR.LT.SNPHYS)SNPHYS=SPAIR IF (SPAIR.LE.STEP)THEN STEP = SPAIR MEC =6 ENDIF ENDIF C C=======> Get next Compton s c a t t e r i n g point C ZINTCO=ZINTCO-(SLENG-SLCOMP)/STEPCO SLCOMP=SLENG STEPC0=GEKRT1*Q(JC0MP+IEKBIN)+GEKRAT*Q(JC0MP+1EKBIN+1) SCOMP=STEPCO*ZINTCO IF(SCOMP.LT.SNPHYS)SNPHYS=SCOMP IF (SCOMP.LE.STEP)THEN STEP = SCOMP MEC = 7 ENDIF C C=======> Get next photo-electron production point C IFCIEKBIN.LT.45)THEN ZINTPH=ZINTPH-CSLENG-SLPHOT)/STEPPH SLPHOT=SLENG STEPPH=GEKRT1*QCJPHOT+IEKBIN)+GEKRAT*QCJPHOT+IEKBIN+1) SPHOT=STEPPH*ZINTPH IFCSPHOT.LT.SNPHYS)SNPHYS=:SPHOT IF CSPHOT.LE.STEP)THEN STEP = SPHOT MEC = 8 ENDIF ENDIF C C=======> Get next p h o t o - f i s s i o n point C IF(JPFIS.GT.O)THEN ZINTPF=ZINTPF-CSLENG-SLPFIS)/STEPPF SLPFIS=SLENG STEPPF=GEKRT1*QCJPFIS+IEKBIN)+GEKRAT*QCJPFIS+IEKBIN+1) SPFIS=STEPPF*ZINTPF IF(SPFIS.LT.SNPHYS)SNPHYS=SPFIS IF (SPFIS.LE.STEP)THEN STEP = SPFIS MEC = 23 ENDIF ENDIF 66 Appendix D. External Conversion code c C C Linear transport 10 20 C C C C C C C c C C 30 IF (EXTCON) STEP=USER STEP EPS = 0.01*EPSIL ~ IFCSTEP.LT.EPS)STEP=EPS DO 20 1=1,3 VECT(I) = VECT(I) + STEP * VECTO+3) SLENG=SLENG+STEP SNPHYS=SNPHYS-STEP DO 30 1=1,7 VOUT(I)=VECT(I) CONTINUE NMEC=1 LMEC(1)=MEC Generate i n t e r a c t i o n IF (MEC.EQ.DTHEN INWVOL = 2 ELSEIF (MEC.EQ.6)THEN CALL GPAIRG ELSEIF (MEC.EQ.7)THEN 1GNEXT = 1 CALL GCOMP ELSE IF(MEC.EQ.8)THEN CALL GPHOT ELSEIF(MEC.EQ.23)THEN CALL GPFIS ENDIF END Appendix E Example Results The following figures show some sample results of a Monte Carlo simulation. The real data is first subjected to geometrical cuts, to eliminate events which do not originate in the target and are most likely external conversions. Multiplicity cuts in the inner two wire chambers also serve to cut clown on external conversions in the target (events with extra low-energy electrons are thus eliminated). A final fairly hard cut in transverse momentum and opening angle ensures that we do not see the detector acceptance (which falls off rather sharply outside these ranges). Figure 11 shows the 7-energy distribution for the real data, calculated by assuming the e e~ pair came from + the decay of a 7r°: E^ E o — Ep { w a r with E„o = 137.8556MeV and Ei pa r being calculated from the 4-momenta of the electrons (information given by the response of the detector). In figure 11, we see that there are two distinct peaks; one centered at E 1 & 0 M e V and the other at E n 50 MeV. The events in the small peak thus correspond to events where no photon was present; TT~ + p —> e + e~ + n + and 7T° -> e + e~ + Appendix E. Example Results 69 while the larger peak contains events in which a photon carried away some energy: the Dalitz decay 0 7r + 1 - 1 —>• e + e +7 with some contamination from the rarer double-Dalitz decay 7T° -> e + e~ + e + e~ + Thus, by making the requirement that E + n should be greater than 20 M e V , we can essentially isolate the Dalitz decay. The simulated data is pure Dalitz decay; no external conversion or other contamination is included. The same cuts are applied as for the real data. Figures 12 and 13 show the comparison of the real data and the simulation results (see equations 2.7 and 2.8 in section 2.1, as well as section 2.3 for a discussion of the variables plotted). Both samples have roughly equal statistics, and have been normalized to the total number of events in each sample. The agreement between data and Monte Carlo simulation is quite good, indicating that the detector response is fairly well modelled, and that contamination of the real data by external conversions and double-Dalitz decays is low. Appendix E. Example Results -20 O 70 20 40 E y 60 80 M e V Figure 11: Real data distribution of the "missing" E-, energy. The ne e~ peak is on the left, the Dalitz peak on the right. The width of the ne e~ peak reflects the resolution of the S I N D R U M detector. The fact that the peak is not centered at 0 is due to imperfect knowledge of the magnetic field. + + Appendix E. Example Results 00 X 71 0.3 (Invariant 0.4 m a s s of e*e" 0.6 pair 0.8 normalized to 10 pion mass) Figure 12: x-distribution for real and simulated data. •olid Una: 19.196 I L_ r«a! Y data I (_ tvtnta (energy partition) Figure 13: y-distribution for real and simulated data.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Monte Carlo simulation of rare decays of the [Pi]⁰
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Monte Carlo simulation of rare decays of the [Pi]⁰ Drees, Reena Meijer 1988
pdf
Page Metadata
Item Metadata
Title | Monte Carlo simulation of rare decays of the [Pi]⁰ |
Creator |
Drees, Reena Meijer |
Publisher | University of British Columbia |
Date Issued | 1988 |
Description | A basic introduction to the techniques of writing computer simulations for particle physics experiments concerning rare decays is given. A very brief outline of the theoretical basis for cross-section calculations is followed by more in-depth discussions of Monte Carlo sampling techniques, including the use of the CERN package DIVON. The reader is introduced to the concepts involved in detector simulation; specifically, the CERN package GEANT in presented and discussed in some detail. Finally, methods of modelling external conversion processes are taken up. Illustrations for the discussions are taken from the various rare decay modes of the π⁰. Example code and some sample results comparing simulated and real data are given in the appendices. |
Subject |
Particles -- Computer simulation |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-08-31 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0084957 |
URI | http://hdl.handle.net/2429/28025 |
Degree |
Master of Science - MSc |
Program |
Physics |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-UBC_1988_A6_7 M44.pdf [ 3.56MB ]
- Metadata
- JSON: 831-1.0084957.json
- JSON-LD: 831-1.0084957-ld.json
- RDF/XML (Pretty): 831-1.0084957-rdf.xml
- RDF/JSON: 831-1.0084957-rdf.json
- Turtle: 831-1.0084957-turtle.txt
- N-Triples: 831-1.0084957-rdf-ntriples.txt
- Original Record: 831-1.0084957-source.json
- Full Text
- 831-1.0084957-fulltext.txt
- Citation
- 831-1.0084957.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
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"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0084957/manifest