Notes on the three-generation neutrino oscillation routines
         -----------------------------------------------------------
                                                             David. A. Petyt
                                                               25-Mar-1996
     The four routines: OSCINIT.FOR, OSCNU.FOR, PROBNU.FOR and PROBNUBAR.FOR
are designed to take the output of a neutrino event generator (neutrino type,
energy and path length) and calculate transition probabilities to all
possible final states, given a mixing scheme specified by the user.
     The structure of the code is as follows:
         1) The user calls the routine OSCINIT.FOR with a set of mixing
            parameters. This is done ONCE, at the beginning of the run.
            *  General three flavour mixing is parameterized by three mixing
               angles and a CP violating phase. OSCINIT uses the standard
               parameterization of Maiani (Phys Lett B62 (1976) p183).
               User supplies three mixing angles    - THETA(1:3).
                             one CP violating phase - DCP
                             (angles are in RADIANS)
               User supplies two delta m^2's - DMSQ(1:2) - units are [eV]^2
            *  Four standard oscillation modes can be selected by the character
               string MODE:
                   MODE                      Oscillation Mode
                   ----                      ----------------
                   `HPS'           3 generation maximal mixing a la HPS
                   `NMT'           2-gen nu_mu->nu_tau with Sin^2(2theta)=1
                   `NME'           2-gen nu_mu->nu_e with Sin^2(2theta)=1
                   `NET'           2-gen nu_e->nu_tau with Sin^2(2theta)=1
               delta m^2 is taken as DMSQ(1) and THETA(1:3),DCP are ignored
               In these cases, the mixing matrix is set up automatically.
            *  There is an option to print out diagnostic information, mixing
               parameters, matrix elements etc. This is controlled by the
               flag IPRINT. If IPRINT=1, then the relevant information is
               written out (useful for MC runs), otherwise it is not.
            *  For more details on how to specify mixing modes in three
               generations, see NuMI-L-129.
         2) The routine OSCNU.FOR is run on EVERY generated neutrino event.
            It returns oscillation probabilities for all possible neutrino
            flavour transitions.
            *  The user supplies :
                  Initial neutrino flavour, FLAVIN, currently defined
                  according to SOUDAN particle code convention:
                         Flavour           Particle code
                         -------           -------------
                          nu_e                    5
                          nu_e bar                6
                          nu_mu                   7
                          nu_mu bar               8
                          nu_tau                  9
                          nu_tau bar             10
                   This can be changed to any convention (e.g GEANT codes)
                   by editing the CODE definitions.
                  Neutrino energy, ENU, in GeV.
                  Neutrino path length, L, in kilometres.
            *  In addition, the user creates and zeroes two 1-dim arrays
               FLAVOUT(6) and OSCPROB(6).
               FLAVOUT and OSCPROB are complementary. FLAVOUT is a list of
               the Soudan particle codes of the 6 possible final state
               neutrino flavours. It is filled in the routine OSCNU.FOR.
               OSCPROB(n) is the transition probability between initial state
               FLAVIN and final state FLAVOUT(n).
            *  OSCNU.FOR determines whether FLAVIN is neutrino or
               anti-neutrino and calls the relevant probability calculating
               function: PROBNU.FOR for neutrinos or PROBNUBAR.FOR for anti-
               neutrinos. All transition probabilities between neutrinos and
               anti-neutrinos are set to zero.
********************************************************************************
Appendix: Obtaining two-generation oscillations with any Sin^2(2theta)
----------------------------------------------------------------------
      The program OSCINIT.FOR provides for the three 2-generation oscillation
modes with Sin^2(2theta)=1. There are occasions, however, when one would
like to simulate two generation oscillations with a smaller value of
Sin^2(2theta) - for example the `Sciulli Point' (Sin^2(2theta)=0.7,
delta m^2=0.01 eV^2). To do this, one has to revert to the three angles in
THETA. Note that for two generation oscillations, there is no CP violation
so DCP is always zero.
      Oscillations in two generations are characterised by one angle, theta,
which is trivially related to the oscillation amplitude Sin^2(2theta). To
obtain a two generation model from the general 3-flavour framework, one
must place theta in one of the elements THETA(n) and set the other two
elements of THETA to zero. One must also place the value of delta m^2 in one
of the DMSQ(n) and set the other DMSQ to zero. The full prescription is as
follows:
--------------------------------------------------------------------------------
|     Mode            | THETA(1)  THETA(2)  THETA(3)  DCP  |  DMSQ(1) | DMSQ(2)|
--------------------------------------------------------------------------------
|                     |                                    |          |        |
|2-gen nu_mu->nu_e    |  theta      0.0       0.0     0.0  |   D m^2  |  0.0   |
|                     |                                    |          |        |
|2-gen nu_e->nu_tau   |   0.0      theta      0.0     0.0  |    0.0   | D m^2  |
|                     |                                    |          |        |
|2-gen nu_mu->nu_tau  |   0.0       0.0      theta    0.0  |    0.0   | D m^2  |
|                     |                                    |          |        |
|3-gen maximal mixing |  0.785     0.615     0.785   1.57  |   Dm21   |  Dm32  |
|                     |                                    |          |        |
--------------------------------------------------------------------------------
* All angles are in RADIANS. Note the placement of the delta m^2's.
* Note that 2-gen nu_mu->nu_e is equivalent to 2-gen nu_e->nu_mu.
* Three generation maximal mixing entered in this way as opposed to MODE='HPS'
  allows one to enter two delta m^2's - Dm21 and Dm32. If both of these are
  of order 10^(-2) eV^2 then potentially large CP violating effects can be
  observed in MINOS (see NuMI-L-129)
*** For example, if one wishes to generate nu_mu->nu_tau oscillations with
the Sciulli parameters (Sin^2(2theta)=0.7,delta m^2=0.01), then the following
values are entered:
                             THETA(1) = 0.0
                             THETA(2) = 0.0
                             THETA(3) = 0.496 radians
                             DCP      = 0.0
                             DMSQ(1)  = 0.0
                             DMSQ(2)  = 0.01
     - where 0.496 is obtained from
   Sin^2(2theta)=0.7 --> Sin(2theta)=0.837 --> 2theta=0.991 --> theta=0.496
********************************************************************************
*                                                                              *
*       Errors/Comments/Questions to  -   d.petyt1@physics.oxford.ac.uk        *
*                                                                              *
********************************************************************************