X-Sun-Data-Type: default
X-Sun-Data-Description: default
X-Sun-Data-Name: simple_3BROC.C
X-Sun-Charset: us-ascii
X-Sun-Content-Lines: 92

/* Simple model of false alarm probability for 3-arm PNS from TSA notes of 8 Aug 97.

   Parameters:   d_psf - diameter of psf in pixels
                 d_lambda - length of lambda search range in pixels
                 npix - total number of pixels in CCD
		 b - background level in electrons per pixel
		 S - signal level in electrons
		 snr - S / sqrt(B + S)
		 thresh - sigma above background to record detection in either PN box or
		      undispersed box
		 f: fraction of incoming light split off into third beam
		 spliteff: efficiency of the beam splitter

		 NOTE THAT B and S are levels that would result for two beam case - the
                      third beam is parasitic off these levels

   Calculation:
                 B = b * d_psf^2
                 B2 = B * (1 - f)
		 B3 = B * f / (0.5 * spliteff)

                 p_fa = npix/d_psf^2 * p_above(B2 + sqrt(B2)*thresh, B2)^2 * d_lambda/d_psf * p_above(B3 + sqrt(B3)*thresh, B3)

		 p_d = p_above(B2 + sqrt(B2)*thresh, B2 + S2)^2 *  p_above(B3 + sqrt(B3)*thresh, B3 + S3)

		 where p_above(n, r) is the probability of drawing n or more counts
		       from a Poisson distribution with expectation r
		      
 */
                 
#include <iostream.h>
#include <stdlib.h>
#include <math.h>
#include <sunmath.h>

extern "C" {double gammq(double a, double x); };

void usage() 
{
  cerr << "simple_ROC <B> <snr> <nsteps>\n";
}

double p_above(double k, double x)
{
  //  cerr << aint(k) << " " << x << "\n";

  return( 1.0 - gammq(aint(k), x));
}

main (int argc, char **argv) 
{
  if (argc != 4) {
    usage();
    exit(1);
  }

  double d_psf = 0.8*2.565;              // 0.8 arcsec seeing, 2.565 pixel/arcsec
  double d_lambda = 40.0;                // +/- 600km/sec at OIII, 2.03 pixel/angstrom
  double npix = 2048*2048;
  double f = 0.259;
  double spliteff = 0.7;

  double b = atof(argv[1]);
  double B = b * d_psf * d_psf;
  double snr = atof(argv[2]);
  double snr_sq = snr*snr;
  double S = (snr_sq + sqrt(snr_sq*snr_sq + 4.0*snr_sq*B))/2.0;

  double B2 = B*(1.0 - f);
  double B3 = B * f / (0.5 * spliteff);
  double S2 = S*(1.0 - f);
  double S3 = S * f / (0.5 * spliteff);
  double snr2 = S2/sqrt(S2 + B2);
  double snr3 = S3/sqrt(S3 + B3);

  cout << "f, snr2, snr3: " << f << " " << snr2 << " " << snr3 << "\n";

  int nsteps = atoi(argv[3]);

  double thresh_min = 0.1;
  double thresh_max = 5.0;
  double thresh_step = (thresh_max - thresh_min)/nsteps;

  double p_fa, p_d;
  for (double thresh = thresh_min; thresh <= thresh_max; thresh += thresh_step) {
    p_fa = npix/pow(d_psf,2.0) * pow(p_above(B2 + sqrt(B2)*thresh, B2),2.0) * d_lambda/d_psf * p_above(B3 + sqrt(B3)*thresh, B3);

    p_d = pow(p_above(B2 + sqrt(B2)*thresh, B2 + S2),2.0) * p_above(B3 + sqrt(B3)*thresh, B3 + S3);

    cout << "B: " << B << "\tsnr: " << snr << "\tS: " << S << "\tthresh: " << thresh << "\tp_fa: " << p_fa << "\tp_d: " << p_d << "\n";
  }
}
----------
X-Sun-Data-Type: default
X-Sun-Data-Description: default
X-Sun-Data-Name: simple_ROC.C
X-Sun-Charset: us-ascii
X-Sun-Content-Lines: 72

/* Simple model of false alarm probability for 2-arm PNS from TSA notes of 8 Aug 97.

   Parameters:   d_psf - diameter of psf in pixels
                 d_lambda - length of lambda search range in pixels
                 npix - total number of pixels in CCD
		 b - background level in electrons per pixel
		 S - signal level in electrons
		 snr - S / sqrt(B + S)
		 thresh - sigma above background to record detection in either PN box

   Calculation: 
                 B = b * d_psf^2
                 p_fa = npix/d_psf^2 * p_above(B + sqrt(B)*thresh, B)^2 * d_lambda/d_psf
		 p_d = p_above(B + sqrt(B)*thresh, B + S)^2

		 where p_above(n, r) is the probability of drawing n or more counts
		       from a Poisson distribution with expectation r
		      
 */
                 
#include <iostream.h>
#include <stdlib.h>
#include <math.h>
#include <sunmath.h>

extern "C" {double gammq(double a, double x); };

void usage() 
{
  cerr << "simple_ROC <B> <snr> <nsteps>\n";
}

double p_above(double k, double x)
{
  //  cerr << aint(k) << " " << x << "\n";

  return( 1.0 - gammq(aint(k), x));
}

main (int argc, char **argv) 
{
  if (argc != 4) {
    usage();
    exit(1);
  }

  double d_psf = 0.8*2.565;              // 0.8 arcsec seeing, 2.565 pixel/arcsec
  double d_lambda = 40.0;                // +/- 600km/sec at OIII, 2.03 pixel/angstrom
  double npix = 2048*2048;

  double b = atof(argv[1]);
  double B = b * d_psf * d_psf;

  double snr = atof(argv[2]);
  double snr2 = snr*snr;
  double S = (snr2 + sqrt(snr2*snr2 + 4.0*snr2*B))/2.0;

  int nsteps = atoi(argv[3]);

  double thresh_min = 0.1;
  double thresh_max = 5.0;
  double thresh_step = (thresh_max - thresh_min)/nsteps;

  double p_fa, p_d;
  for (double thresh = thresh_min; thresh <= thresh_max; thresh += thresh_step) {
    p_fa = npix/pow(d_psf,2.0) * pow(p_above(B + sqrt(B)*thresh, B),2.0) * d_lambda/d_psf;

    p_d = pow(p_above(B + sqrt(B)*thresh, B + S),2.0);

    cout << "B: " << B << "\tsnr: " << snr << "\tS: " << S << "\tthresh: " << thresh << "\tp_fa: " << p_fa << "\tp_d: " << p_d << "\n";
  }
}