JNs.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <cmath>
 
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
 
#include "JLang/JManip.hh"
#include "JMath/JMathSupportkit.hh"
 
#include "Jeep/JParser.hh"
 
 
/**
 * Auxiliary program to determine number of signal events for given background and number of standard deviations.
 */
int main(int argc, char* argv[])
{
  using namespace std;
  using namespace JPP;
 
  string    outputFile;
  double    stdev;
  double    precision;
  int       debug;
 
  try {
 
    JParser<> zap("Auxiliary program to determine number of signal events for given background and number of standard deviations.");
 
    zap['o'] = make_field(outputFile)          =  "poisson.root";
    zap['s'] = make_field(stdev);
    zap['p'] = make_field(precision);
    zap['d'] = make_field(debug)               = 2;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  const double STDEV_MIN = 1.0;
  const double STDEV_MAX = 7.0;
 
  if (stdev < STDEV_MIN || stdev > STDEV_MAX) {
    FATAL("Number of standard deviations out of range [" << STDEV_MIN << "," << STDEV_MAX << "]" << endl);
  }
  
  const double P = getP(stdev);
 
  cout << "probability (single-sided Gauss): " << SCIENTIFIC(12,3) << P << endl;
 
 
  TFile out(outputFile.c_str(), "recreate");
 
  TH1D h0("h0", NULL, 500, -2.0, 2.0);
  TH1D h1("h1", NULL, 500, -2.0, 2.0);
 
  for (Int_t ix = 1; ix <= h1.GetXaxis()->GetNbins(); ++ix) {
 
    const double x = pow(10.0, h1.GetXaxis()->GetBinCenter(ix));
 
    const size_t n = getNs(x, P);
    const double f = getFs(x, P, precision);
    
    DEBUG("    x = " << FIXED(7,3) << x
          << " n = " << FIXED(7,3) << n - x 
          << " f = " << FIXED(7,3) << f - x << endl);
 
    h0.SetBinContent(ix, n - x);
    h1.SetBinContent(ix, f - x);
  }
 
  out.Write();
  out.Close();
}