M O N T E C A R L O S I M U L A T I O N O F R A R E D E C A Y S O F T H E 7T° Reena Meijer Drees B. Sc. (Hons, co-op) University of Waterloo A T H E S I S S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F M A S T E R O F S C I E N C E in T H E F A C U L T Y O F G R A D U A T E S T U D I E S D E P A R T M E N T O F P H Y S I C S We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A September 1988 (c) Reena Meijer Drees In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department The University of British Columbia Vancouver, Canada DE-6 (2/88) Abs t r ac t 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 the-oretical 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 i i Acknowledgement vi 1 Introduction 1 2 Data 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 External Conversions and Compton Scattering 38 4 Conclusions 43 Bibliography 45 A D I V O N coded example 47 i n B G E A N T detector defini t ion 50 C G E A N T t r ack ing code 56 D E x t e r n a l Convers ion code 60 E E x a m p l e Resul t s 68 i v Lis t 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 mod-elling detector acceptance 21 6 Z distribution of event vertices; dotted line shows result of modelling detector acceptance 21 7 Simplified G E A N T program structure (taken directly from the GEANT User's Manual) 31 8 Total cross section for pair production (barn/atom) as calculated using G E A N T parameterized fits of cross section data 40 9 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 (NSERC) 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 Mr. C. Niebuhr, for their help and suggestions. v i Chap te r 1 In t roduc t ion The collection of the theories of weak interactions, quantum electrodynamics (QED) and quantum chromodynamics (QCD) 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. Sup-pressed 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 non-conventional theories by virtue of their sensitivity to "extra" physics. However, their very rareness makes them difficult to study experimentally. Measurements of rare de-cay 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. Wi th 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 simu-lation 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 exper-imental 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 . Chap te r 2 D a t a S imula t ion : Event Genera t ion 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, x2,..., xn) 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 gener-ation 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 4 2.1 T h e Different ial 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 distribu-tion of the e+e" pair, etc. This information is given in a quantity called the "differential cross section" 1 dcr(xi,X2, • • •, xn), 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,..., xn) dxi • • • dxn — 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 me-chanics 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>n orthonormal): Ho^n = En(f)n 1 or "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 unper-turbed system: To find the coefficients an we multiply through by (j>*j and integrate, using the orthonor-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 af(t) = -i J* dt' j d3x <f>) V fa e ' ( s / -SO*' 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 Tfi = af(T/2) = -if dt [ d3x[4>j(x)e-iE^}*V(x,t)[U^y-lEtt} J-T/2 J = -i J dAx<t>){x)V{x)(j)l{x) 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 poten-tial, we consider Maxwell's equations and make the substitution pA< _» 7y< + eA» 2 Thi 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 Wu = l im Chapter 2. Data Simulation: Event Generation 6 Figure 1: Feynman diagram Figure 2: Feynman diagram for 7 emission for 7 exchange so that the Dirac equation becomes with the perturbation Carrying on as in the Schrodinger case, we obtain Tfi = -iJ jfA^x (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: • 2 A " = ft Chapter 2. Data Simulation: Event Generation 7 with solution Q where q = pa — pc, the momentum carried away by the photon. Figure 2 shows the corresponding Feynman diagram. Substituting the result above into equation 2.2 leads to Tfi = -iJ*xjl (-^j to first order. Performing the x integration, (recall the explicit definition of the tran-sition current j M , equation 2.3) we get TH — (-euc-y^a) ( -eM D 7 M u i )(27r) 4 8(pa + pb - pc - pd) = -i (27r)4 8(pa +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/q2 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 (2.4) r where F and dQ symbolically represent the incident flux and the Lorentz invariant phase space factor; 3 F= \va\2Ea-2Eb where Pa va = —f- = velocity of beam particle 3for 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 ) 4 < ^ + W - Pi - • • • - Pn) (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 an-other force called the "strong force", for which a candidate theory is quantum Chapter 2. Data Simulation: Event Generation 9 chromodymamics (QCD). Since QCD is not as calculable as Q E D , the usual tech-nique 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. - r - c £ * / ; * ^ K i - 3 , ( i + » , + ^ ) (2-6) where m is the electron mass and / i is the mass of the 7r° , and 4?n2 ~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 (2-7) y = energy partition = -rz, zri (2-8) \P++P-\ 2.2 M o n t e C a r l o Sampl ing M e t h o d s The kinematics of each event are chosen randomly according to the differential branch-ing 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~ax, we solve the integral equation u= f e - a x ' d x ' Jo where u is a random number on the interval [0,1), for x, to obtain x = ln( 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 x0 lies in the domain of f(x). F(xQ) is evaluated, and then a new random number y is chosen, where y lies in the range of F. The random number x0 is kept if y < F(x0). Note that this is an exact technique. 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 x2, and so on. We again need the cumulative probabilities, integrating over the quantities we are not choosing: Once Xi 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 Eu-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.4 Extensive useful scientific literature is also available from CERN's Scientific Information Service. 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, 4further information obtainable from: Program Library, Division D D , C E R N , 1211 Geneva 23, Switzerland D I V O N 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 DIVON'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 sub-regions 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,Si) ^ S P R D M A X 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 s function 5 C O M M O N / S A M P L E / N P O I N T 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. In tegrat ion: D I V O N supports several alternative numerical integration proce-dures: • quasi-uniform Monte Carlo estimate (recommended) • pseudo-random estimate with importance sampling • 2 n d , 3 r d , or 5 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. Sampl ing : 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 rec-ommended 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 parti-tioning; 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. Phase 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 L a b o r a t o r y Frame K i n e m a t i c s The random kinematical quantities having been chosen according to the differential branching ratio, the event can now be materialized in the laboratory. From the cho-sen 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 (2.10) + 7 (2-11) + 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 MeV fa = 0.2037 En = 939.9724 MeV (3n = 0.0299 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 px, py and pz randomly, which, expressed mathematically, means sampling the distribution dpx dpy dpz subject to the constraint 7T —> e + e -»• e + + e~ -» e + + e~ Chapter 2. Data Simulation: Event Generation 17 Converting to polar coordinates, this is equivalent to sampling J j d<f> rf(cos 6) 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: Elab = jE -j/3 -p Pub = p-MiE-ft-p-1 with P = 7 1 7 + 1 PIT0 Eir° 1 the usual Lorentz quantities. For reaction 2.13 sampling is performed to obtain the kinematics from the differ-ential branching ratio. In this case, we use the results from Krol l 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: Elt - m2n + x2 ^pair — Q jp 2Etot E+ — ~^(Epair ~\~ Impair ' J / ) E- = Epair — E+ Chapter 2. Data Simulation: Event Generation 18 with the angle between the e + and e given by 2E,E_ - 2m2e - x2 cos a — where Etot = ET- + Ep = mv- + mp (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 9n cos <f>n — sin 0n sin 6n cos <j>n \ TZ — cos 9n sin <f>n cos 8n sin 9n sin <f>n y — sin 6n 0 cos 9n 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 T h e Event Ver t ex To fix the event completely we must localize it inside the detector. This involves gen-erating 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 experi-ment. 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 mate-rial 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 postu-lates 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 data d o i t e d l ine: d i g i t i z a t i o n r e s u l t vertex R (mm) Figure 3: Radial distribution of event vertices I i I so l id l ine , r e a l d a t a d o t t e d l ine: d i g i t i z a t i o n r e s u l t - 1 2 0 Figure 4: Z distribution of event vertices Chapter 2. Data Simulation: Event Generation 21 .025-.020 • solid line: real data dotted line: monte carlo 2f> vertex R (mm) Figure 5: Radial distribution of event vertices; dotted line shows result of modelling detector acceptance. T3 <U £ .016 10 E i— o —> .010 c 3 o s o l i d l ine : r e a l d a t a d o t t e d l ine: m o n t e c a r l o v e r t e x Z ( m m ) 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 low-mass 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, 6 the technique is to force 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. 6i.e. stepping the 7 through the detector, at each step calculating the probability for conversion pc and then allowing conversion if (random number) < pc. Obviously if pc is small we will need many 7 ' s to see an event. 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 ex-ternal conversions go as e, then the probability for a photon to undergo two conversions goes as e2. If one now increases e by a factor of 10 to make conversions more likely, 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 sam-pling 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~ + 7 (1.198 x 10" 2 branch) (2.14) TT° -> e+ + e~ + e+ + e~ (3.24 x 1 0 - 5 branch) (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 conver-sion 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 tepping : 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 el-ements (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 examin-ing 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. • Ou tpu 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. • Graph ic 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 Chapter 3. Detector Simulation: Tracking 26 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 : In t roduc t ion and Genera l Ph i losophy 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 ZE-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 UGINIT, • 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. In i t i a l i za t ion : 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: Tracking 28 © GINIT—initializes G E A N T common blocks with default options (includ-ing 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—init ial izes the graphics package (if desired). G E A N T also needs to initialize its tables of particle properties, material proper-ties, 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 (radia-tion 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 hydro-gen, deuterium, helium and neon are for the liquid state. The appropriate values for gaseous helium.are: density = 1.78 x 10" 4 g cm'3 radiation length = 1.075120 x 106 cm absorption length = 6.80 x 105 cm 1 Note 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 ( IBREM=1) , and if this is not desired, the user should include the data card B R E M 2 in the input. Chapter 3. Detector Simulation: Tracking 29 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 preci-sion. 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 pre-viously 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. Genera t ion and Track ing : 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: Chapter 3. Detector Simulation: Tracking 30 — 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 accumu-late 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. Te rmina t ion : 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 Detec tor Def in i t ion 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 UGINIT 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, positioned in one another. The user defines an initial mother volume, and positions daughter volumes inside this. These Chapter 3. Detector Simulation: Tracking 31 MAIN (user) — GZEBRA initialisation of ZEBRA ivstcm, dynamic core allocation — UGINIT (user) — GINTT initialisation of GEANT3 variables — GFFGO interpretation of data cards — GZINIT initialisation of ZEBRA core divisions and link areas — GPART creation of the 'particle' data structure JPART — G M A T E creation of the 'material' data structure JMATE — 'user code' description of the geometrical setup, of the sensitive detectors creation of data structures JVOLUM, 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 JXYZ - 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: Tracking 32 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 Mas-ter 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 of the mother into several identical volumes with the same parameters as the mother. 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 GSPOS, and di-vided 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 2 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. 3 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 . 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 associ-ated 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 UGINIT, 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 event-tracking 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 GTNINO (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 mag-netic 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. Chapter 3. Detector Simulation: Tracking 36 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 D i g i t i z a t i o n and Ou tpu t G E A N T provides the user with a set of data banks for recording "hits", the (user-defined) information regarding the track parameters (momentum, energy loss, co-ordinates, 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 exper-iment, 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 GSAHIT can be used to store any desired hit information. GSAHIT should be called after every tracking step, i.e. in G U S T E P . GSAHIT 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 GUDIGI , 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 GFHITS to unpack the information stored in the hit bank using GSAHIT. 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 Chapter 3. Detector Simulation: Tracking 38 its number and the energy deposited in it after every step. In GUDIGI 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 Conversions and C o m p t o n Scat ter ing 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. 4Then, the 4except in the case of rate calculations, where one is estimating the total number of events, regardless of topology, that the detector will see. Chapter 3. Detector Simulation: Tracking 39 relative probability pxc for each of these sections to initiate a photon conversion must be 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 ( C 5 . H 4 O 2 ) f°r a photon of energy 50 MeV. In general, the probability p,- of a particle to undergo an interaction i is given by L T A Pi = T = L- — 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 MeV 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. By 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 first-generation 7 's; second-generation conversions produce electrons of too low energy. Chapter 3. Detector Simulation: Tracking 40 -Sl2 ? o £ 8-c o * 6 -B 5 o •024 o .021 - .018 - .015 tr liquid hydrogen - .012 o 3 .009 — h .006 .o c - .003 % 0. 1 50 — I 150 .000 2 100 200 Gamma energy (MeV) 250 300 Figure 8: Total cross section for pair production (barn/atom) as calculated using G E A N T parameterized fits of cross section data . r 80 40 60 Gamma energy (MeV) 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: Tracking 41 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 STEP, 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. Chap te r 4 Conclus ions 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 exper-imental 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 simu-lation. 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. Bib l iog raphy [1] C. Cohen-Tannoudji, B. Diu, 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 . Krol 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] DIVON4 (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, DD/EE/84 -1 (1986) 45 Bibliography 46 [15] Particle Data Group, Review of Particle Properties: Phys. Let. 170B (1986) and also the associated Particle Properties Data Booklet A p p e n d i x 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, XVECT0R(2), XMIN(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) GOTO 10 C i n i t i a l i z e DIVON i n t e g r a t i o n p a r a m e t e r s , l i m i t s on X XMIN(1) = XMINIMUM XMAX(l) = XMAXIMUM C l i m i t s on Y XMIN(2) = -SQRT(1. - RR) XMAX(2) - -XMIN(2) C c o n t r o l parameters OPEN (UNIT=10, READONLY, SHARED, STATUS='OLD') READ ( 1 0 . 101) SPRDMAX, MAXPTS 101 FORMAT(F10.0, 110) C p a r t i t i o n i n t e g r a t i o n space CALL PARTN(N, XMIN, XMAX, SPRDMAX, MAXPTS) C e s t i m a t e i n t e g r a l READ (10 , 102) JDEG, NPT 102 FORMAT(2I10) CALL INTGRL(N, JDEG, NPT, ANS, ERROR) CLOSE(UNIT=10) INIT = .TRUE. C f i l l a r r a y XVECTOR ( v a l u e s of X, Y) with random X,Y 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, XVECTOR) RETURN END Appendix A. DIVON coded example DOUBLE PRECISION FUNCTION DFUN(N, XY) c*********************************************************** C joint probabi l i ty fn . for X and Y, from Krol 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 / PI, RR, GP, GN, ASLOPE DATA CONST / 5.80857D2 / X = XY(1) Y = XY(2) IF (X . L T . RR) GOTO 100 ! redundant but safe IF (X . L T . XMIN .OR. X .GT. XMAX) GOTO 100 ETA = SQRTC1. - RR/X) IF (ABS(Y) .GT. ETA) GOTO 100 EPAIR = MPO/2. * (1. + X) PPAIR = MPO/2. * (1. - X) EP0 = (EPAIR + PPAIR*Y) / 2. EE0 = EPAIR - EP0 IF (EP0 . L T . EMIN .OR. EE0 . L T . EMIN) GOTO 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 100 DFUN = 0.0D0 RETURN END A p p e n d i x B G E A N T detector defini t ion 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 SCINT is used to define the tracking media POLYSTYRENE and P0LYH0D0. P0LYH0D0 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 51 Figure 10: The S I N D R U M detector. t 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 ASCI NT/12.011,1.008/AKAPA/12.,1.,14.,16./ DATA Z S C I N T / 6 . , 1 . / Z K A P A / 6 . . 1 . . 7 . , 8 . / DATA USCINT/1. ,1.1/WKAPA/19. ,14. ,2. ,4 . / DATA AMYLAR/12..1. .16./AMAKRO/12.,1. ,16./ DATA ZMYLAR/6 . ,1 . ,8 . /ZMAKR0/6 . ,1 . ,8 . / DATA WMYLAR/5.,4.,2./WMAKRO/3..8.,2./ DATA CHANAM/'CHA1','CHA2','CHA3','CHA4','CHA5','CHA6','CHA7'/ DATA GASNAM/'ARG1','ARG2','ARG3','ARG4','ARG5','ARG6','ARG7'/ DATA ACTNAM/'ACT1' , 'ACT2' , 'ACT3' , 'ACT4' , 'ACT5' , 'ACT6' , 'ACT7' / DATA UIRNAM/'WIR1','wIR2','UIR3','WIR4','WIR5','wIR6','wIR7'/ C read in SINDRUM geometry from EGEO bank CALL READ_EGEO C Define materials 1.01, 1.,0.071 .865..790.,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) 4.00, 2.,0.00000178,1.E8,1.E8,0,0) CALL GSMATE( 1,'HYDROGENS CALL GSMATE( 3,'HELIUMS CALL GSMATEC 9,'ALUMINIUMS CALL GSMATE(10,'I RONS CALL GSMATEC11,'COPPERS CALL GSMATE(15,'AIR$ C f ix up GEANT bug with vacuum CALL GSMATEC16,'VACUUMS ' , CALL GSMATE(17,'ARGONDENSS', 39.948,18.,0.00178,16500.,0.,0,0) CALL GSMIXTC20,'SCINTS ',ASCI NT,ZSCINT,1.032,-2,WSCI NT) CALL GSMIXT(21,'KAPACELLS ',AKAPA,ZKAPA,0.15,-4.WKAPA) CALL GSMIXT(22,'MYLARS ',AMYLAR,ZMYLAR,1.39,-3,WMYLAR) CALL GSMIXT(23,'MAKROLONS ',AMAKRO,ZMAKRO,1.18,-3,WMAKRO) C Defines USER tracking media parameters FIELDM = 3.313 !FIELD STRENGTH IN kG I FIELD = 3 !FIELD CONSTANT ALONG AXIS 3 TMAXFD = 10.0 !MAX. ANGLE DUE TO FIELD (DEGREES) STMIN = 0.0000001 !MIN. STEP DUE TO MULT. SCATTERING (CM) DMAXMS = 0.5 !MAX. DISPLACEMENT DUE TO MULT. SCATTERING (CM) DEEMAX = 0.2 !MAX. FRACTIONAL ENERGY LOSS EPSIL = 0.01 !TRACKING PRECISION (CM) CALL GSTMED( 1,'DEFAULT MEDIUM HELIUMS',3,0,1 FIELD,FIELDM, * TMAXFD,DMAXMS,DEEMAX,EPSIL,STMIN, 0,0) DMAXMS = 0.01 !MAX. DISPLACEMENT DUE TO MULT. SCATTERING (CM) DEEMAX = 0.003 !MAX. FRACTIONAL ENERGY LOSS EPSIL = 0.0003 ITRACKING PRECISION (CM) Appendix B. GEANT detector definition 53 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 = 0.5 !MAX. DISPLACEMENT DUE TO MULT. SCATTERING CCM) DEEMAX = 0.2 !MAX. FRACTIONAL ENERGY LOSS EPSIL = 0.01 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 cyl inder: PARC1)=35.6 PARC2)=38.4 PAR(3)=66.0 CALL GSV0LUC'ALUM','TUBE',A,PAR,3,IV0L ) C magnet co i I : PARC1)=38.4 PARC2)=40.7 PARC3)=66.0 CALL GSV0LUC'C0IL','TUBE',5,PAR,3,IV0L ) Appendix B. GEANT detector definition C Posit ion 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 GSPOS('HODO',1, 'SIND',0. ,0. ,0. ,0, 'ONLY') CALL GSPOS('ALUM',1, 'SIND',0 . ,0 . ,0 . ,0 , 'ONLY') CALL GSPOS('COIL' ,1 , 'SIND' ,0 . ,0 . ,0 . ,0 , 'ONLY') 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, (obligatory 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 rotations for various runs REAL DPHINEU(5, 8) INTEGER NEUGEO(8) DATA NEWGEO / 0, 92, 160, 171, 208, 225, 265, 270 / DATA DPHINEW / 1 0.41425, 1.69182, 0.1199. 0.28887. 0.17578. 2 0.45825, 1.72272, -0.01852, 0.22555, 0.17578, 3 -0.259 , 1.71401, 0.04749, 0.09988, 0.17578, 4 2.52689, 1.71401, 0.04749, 0.09988, 0.17578, 5 1.404 , 1.68237, 0.02747, 0.0992, 0.17578. 6 1.9934 , 1.71338, 0.04483, 0.09966, 0.17578, 7 2.21378, 1.97664, 0.20499, 0.39048, 0.17578, 8 2.26254, 2.00748, 0.46625, 0.53028, 0.17578 / C f ind out geometry for this run DO 111=1,8 IF(ABS(NRUN).LE.NEWGEO(IU)) GOTO 10 END DO 10 IU=IU-1 WRITE(6, 1000) IU 1000 FORMATC ROTATION SET USED:', 12 , / , / ) Appendix B. GEANT detector definition C f i l l geometry commons with data from bank EGEO (for C chambers and hodoscope) C ATTENTION: only one hodoscope supported; angle between C anode and inner cathode is assumed to be zero as well C as z-posit ion of middle of chamber for C 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 ! HODOSCOPE GEOMETRY NSECT=IDATA(IP+2) HODPHI0=RDATA(IP+3)*RADDEG HODLEN=RDATA(IP+4)/10 HODTHI=RDATA(IP+5)/10 IP=IP+LTRHOD 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 !OUTER CATHODES 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 A p p e n d i x C G E A N T t rack ing 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 SINDRUM), 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 Called 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 !number of tracking s NUMTRA=NTRACK DO I=1,NUHTRA NHITS(I)=0 HODOREACHED(I)=.FALSE. ENDDO CALL BDROP(IDATA,'E') !clear out hit bank CALL BGARBCIDATA) CALL GTREVE !do the tracking RETURN END Appendix C. GEANT tracking code SUBROUTINE GUSTEP C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C User routine ca l led at the end of each tracking step C INUVOL is different from 0 when the track has reached C a volume boundary C ISTOP is different 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 NACTUAL=NACTUAL+1 IFCNACTUAL .GT. NSTEPS)THEN WRITE(6,900) IDEVT 900 FORMAT(2X, 15, ' EVENT ABORTED - TOO MANY TRACKING STEPS') IEOTRI=1 RETURN END IF C store secondaries into stack - - only those originating C inside target or chambers IF(NGKINE.GT.O) THEN IF(TEXTCNLEVEL) .NE. 'COIL' .AND. TEXT(NLEVEL) .NE. '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 sensit ive volume—chambers or hodos only? IF(.NOT.(INCHAMBER.OR.INHODOSKO))RETURN C 100 hits only IF(NHITS(ITRA).GE.100)RETURN C only interested in e+ or e- part ic les in gas IF(INCHAMBER.AND.ITRTYP.NE.2)RETURN C just entering volume? IF(INUVOL.EQ.1)THEN C store number of hi ts in this 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 exit ing volume? IF(INUVOL.GE.2.OR.ISTOP.NE.0)THEN C then store Pz at exit point, 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 A p p e n d i x D E x t e r n a l Convers ion 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 7T° -> 7 + 7 (D.16) 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 user-provided USER_STEP. 60 Appendix D. External Conversion code Q*********************************************************** SUBROUTINE GUTRAK C User routine to control tracking at track leve l . 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 is this 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 (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 = 3 d 3 H n i ( i 1 f 1 i , - 0 3 ' ( £ : ) ( 1 3 A 3 1 N)iX31 ) d I I = 3 a 3 H n I ( , H V l , - 0 3 - ( £ : ) O 3 A 3 1 N)lX31 ) d I 0 = SaSHt t l 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 • 3 n b l ' = N031X3 1NV33 JOJ- 6e i j . uo.isjaAUOo i s u j a j x a }as os '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 0 3(r3N - d A i H i r a 0 T 3 N - 1 0 A r l N I ) d I 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 seq JUSAG am >(33qo 3 dI0N3 N a n i s a 1 = i a i o s i (,Sd31S ONIXOVMl ANVW 001 - 0 3 i a o a V < 'X£'9I'/1N3A3 /HVWHOd 006 N3Hi(Sd31SN "13" "IVni3VN)dI l+1Vni3VN = 1VP13VN sd3JS AUBUI 001 3 j . i JU3A3 u o q s - s d s j s 6 u i 5 | 3 e j j j.o jsquinu junco 3 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 / a o d - s i i H 3 : a i a $ 3 a n i 3 N i ; a a m o N i /aod -wvH33:a ia$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 * psddojs seq yoeji j.t 0 WOJJ. j u a j a j j i p s j dOlSI 3 » Ajepunoq a d i n i oA e 3 » P3U0B3J SBU, >|3BJ1 3U,} U3L|M 0 U10J}. 1U3J3J.JIP SI 10AMNI 3 » d s j s 6u_i>|3BJi ipsa j.o pua ai|i IB p a n s a auui ioj j s s n 3 ***********************************************************3 d3 i sn3 B N i i n o a s n s ***********************************************************3 29 9poo UOISI9AUOQ jvujQixg (j xxpusddy 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 interaction length, code plagiarized 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 calculate step size and interaction probabi l i ty 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 else: 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 interaction to be 1.0E-7 * interaction 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 originating C inside target or chambers 1 IF(NGKINE.GT.O) THEN IF(TEXTCNLEVEL) .NE. 'COIL' .AND. TEXT(NLEVEL) .NE. 'ALUM' .AND. TEXT(NLEVEL) .NE. 'SECT') CALL GSKING(O) ENDIF C debug requests IF(ISWIT(3).EQ.DTHEN CALL GSXYZ CALL GPCXYZ END IF C produce hit 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 IF(NHITS(ITRA).GE.100)RETURN IF(INCHAMBER.AND.ITRTYP.NE.2)RETURN IF(INWVOL.EQ.DTHEN NHITS(ITRA) = NHITS(ITRA)+1 NH = NHITS(ITRA) ZENTRY(ITRA,NH) = VECT(3) NUMWIR(ITRA,NH) = NUMBER(NLEVEL)-1 EDEPO(ITRA,NH) = 0. ! HODOSCOPE)? 1ALLOU ONLY 100 HITS !IGNORE PART. EXCEPT e+-!ENTERING VOLUME IWIRES START AT 0 IN GAS 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 £*********************************************************** SUBROUTINE GTGAMA C. C ****************************************************** c! * * C. * A photon track is extrapolated by one step * C. * * C. * ==>Called by : GTVOL * C. * Authors R.Brun, M.Maire ********* * C. * * C. * 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 KCASE,NGKINE REAL 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 Is the part ic l e below threshold energy ? C IF(GEKIN.LE.CUTGAM)THEN ISTOP=2 DESTEP=DESTEP+GEKIN Appendix D. External Conversion code 66 STEP=EPSIL SLENG=SLENG+EPSIL NMEC=1 LMEC(1)=30 RETURN ENDIF C C Update local pointers i f medium change C 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 pair 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 scattering 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 photo-fission 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 Appendix D. External Conversion code c C Linear transport C 10 IF (EXTCON) STEP=USER STEP EPS = 0.01*EPSIL ~ IFCSTEP.LT.EPS)STEP=EPS DO 20 1=1,3 20 VECT(I) = VECT(I) + STEP * VECTO+3) SLENG=SLENG+STEP SNPHYS=SNPHYS-STEP DO 30 1=1,7 VOUT(I)=VECT(I) 30 CONTINUE C NMEC=1 LMEC(1)=MEC C C Generate interaction C C C C c C C 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 A p p e n d i x E E x a m p l e Resul ts 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^ Ewo — Epa{r with E„o = 137.8556MeV and Epair 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 E1 & 0 MeV and the other at En 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 + 1 - 1 7r —>• e + e + 7 with some contamination from the rarer double-Dalitz decay 7T° -> e + + e~ + e + + e~ Thus, by making the requirement that En should be greater than 20 MeV, we can essentially isolate the Dalitz decay. The simulated data is pure Dalitz decay; no external conversion or other contamina-tion 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 70 - 2 0 O 20 40 60 80 Ey 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 71 0 0 0.3 0.4 0.6 0.8 1 0 X ( I n v a r i a n t m a s s o f e*e " p a i r n o r m a l i z e d t o p i o n m a s s ) Figure 12: x-distribution for real and simulated data. I L_ I (_ • o l i d U n a : 19.196 r « a ! d a t a t v t n t a Y ( e n e r g y p a r t i t i o n ) 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 |
AggregatedSourceRepository | 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:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0084957/manifest