/* loudest.c
 * Jolien Creighton <jolien@weyl.phys.uwm.edu>
 * 23 July 1999
 *
 *
 * Compilation and Running:
 *
 * Requires Numerical Recipes C routines ran1() and poidev().  To compile:
 *
 * gcc loudest.c -o loudest -L/usr/local/recipes/lib -lrecipes_c -lm
 *
 * The program takes no arguments.
 *
 *
 * Synopsis:
 *
 * This program performs a Monte Carlo calculation of the probability that
 * the rate of events exceeds the rate limit of Allen et al. or Finn given the
 * efficiency of the loudest event observed.  (Changing between the two rate
 * limits requires modifying the code and recompiling.)
 *
 *
 * Strategy:
 *
 * 1. Choose a population event rate mu from a distribution that is uniform
 *    between zero and mumax.
 *
 * 2. Using this rate, find the number of events m that occur.
 *
 * 3. For each of the m events that occurred, choose a random efficiency eps.
 *
 *    The efficiency of an event is uniformly distributed between zero and one
 *    and represents the probability of another event being louder than this
 *    particular event.  (Loud events have small efficiencies.)
 *
 * 4. Find the minimum efficiency epsmin, i.e., the efficiency of the loudest
 *    of the m events.
 *
 * 5. If the minimum efficiency epsmin is not within some slop amount deps
 *    of the observed efficiency epsobs, then discard this trial (because
 *    it does not agree with the observations) and go back to step 1.
 *
 * 6. Otherwise, increment the total number of acceptable trials ntot,
 *    and, if the rate mu is less than the rate limit mu90, which is obtained
 *    from the formula of Allen et al. or Finn, increment the number of
 *    acceptable trials nyes that satisfy the rate limit as well.
 *
 * 7. Repeat until there are enough acceptable trials.
 *
 * The fraction frac = nyes/ntot should then be 0.9 if the rate limit is
 * correct.  Note that the case when there are zero events in a trial is
 * automatically excluded because the efficiency will not equal the observed
 * efficiency in this case.
 *
 * Note: frac is what I consider to be P(mu > mu90 | epsobs).  Please let me
 * know if this is not so!
 *
 *
 * Results:
 *
 * For mu90 given by Allen et al.:
 *
 *    ntot    nyes    frac     err
 *  262144  235998   90.03%   0.06%
 *
 * For mu90 given by Finn:
 *
 *    ntot    nyes    frac     err
 *  262144  218741   83.44%   0.07%
 *
 * This seems to indicate that the confidence level Finn's upper limit
 * is less than 90%.
 *
 *
 * Systematic Errors:
 *
 * Ideally, we want mumax -> infinity and deps -> 0.  However, the number
 * of trials rejected increases for large mumax and small deps.  The values
 * I have chosen here will bias the result, but hopefully not by much.  I have
 * not tried to quantify these systematic errors.
 * 
 * 
 * References:
 *
 * Allen et al.(1999), Phys Rev Lett (to appear)
 * Finn (1999), gr-qc/9907070
 *
 */

#include <stdio.h>
#include <math.h>

float ran1(long *);
float poidev(float, long *);

int main()
{
  const int   nmax   = 262144; /* maximum number of trials          */
  const float mumax  = 50;     /* maximum rate (should be infinity) */
  const float epsobs = 0.33;   /* observed loudest event efficiency */
  const float deps   = 0.01;   /* allowed slop in efficiency        */

  float mu90 = 3.890/epsobs;   /* rate limit of Allen et al.        */
  /* float mu90 = 9.823; */    /* rate limit of Finn                */
  float frac;                  /* fraction of trials with mu < mu90 */
  long  nyes = 0;              /* number of trials with mu < mu90   */
  long  ntot = 0;              /* total number of trials            */
  long  seed = -101;

  printf("  ntot\t  nyes\t  frac\t   err\n");

  while (ntot < nmax) {
    
    float epsmin = 1;                      /* efficiency of the loudest event */
    float mu     = mumax*ran1(&seed);      /* population rate for this trial  */
    int   m      = (int)poidev(mu, &seed); /* number of events in this trial  */
    
    /* compute the efficiency for the loudest event
     * (the minimum of the efficiencies of all the events that occur)
     */
    while (m-- > 0) {
      float eps = ran1(&seed); /* efficiency of the present event */
      if (eps < epsmin)
        epsmin = eps;
    }
    
    if (fabs(epsmin - epsobs) < deps) {
      /* this trial is a candidate (i.e., it reproduces the observed
       * loudest event efficiency to within the allowed slop)
       */
      ++ntot;
      if (mu < mu90)
        ++nyes;
      if (ntot%64 == 0) {
        frac = (float)nyes/ntot;
        fprintf(stderr, "%6ld\t%6ld\t%6.2f%%\t%6.2f%%\r",
                ntot, nyes, 100*frac, 100*sqrt(frac*(1 - frac)/ntot));
      }
    }
    
  }

  frac = (float)nyes/ntot;
  printf("%6ld\t%6ld\t%6.2f%%\t%6.2f%%\n",
         ntot, nyes, 100*frac, 100*sqrt(frac*(1 - frac)/ntot));

  return 0;
}