// File:  runSigCalc.cc
// Glen Cowan
// RHUL Physics
// 14 April 2008

// Simple driver program for the function getSignificance, which 
// returns the significance with which one can reject a hypothesized
// signal strength parameter mu.  (For discovery, one rejects mu = 0.)

// Input:  runSigCalc reads in a text file with the input data.  
// Lines beginning with # are comments.
// The first (non-comment) line contains the target luminosity.
// The second (non-comment) line contains the accepted number of
// signal events from MC and the luminosity of the signal MC sample.
//
// On subsequent lines the user lists for each background component the
// accepted number of background events m and the effective luminosity 
// L_MC used to generate the background sample.

#include <iostream>
#include <fstream>
#include <cstdlib>
#include <string>
#include <sstream>
#include <vector>
#include "SigCalc.h"
#include "getSignificance.h"

using namespace std;

int main(int argc, char **argv) {

// Read input file

  string inputFileName;
  if ( argc == 1 ) {
    cout << "Enter name of input file: ";
    cin >> inputFileName;
    cin.ignore();         // clear buffer
  }
  else {
    inputFileName = argv[1];
  }

  ifstream f;
  f.open( inputFileName.c_str() );
  if ( f.fail() ){
    cout << "Sorry, couldn't open input file" << endl;
    exit(1);
  }

  double ns, s, lumiTarget, lumiSig;
  vector<double> m;
  vector<double> b;
  vector<double> tau;
  vector<double> w;
  istringstream instream;       // Declare an input string stream
  string line;
  int lineNum = 0;
  while ( getline(f, line) ) {
    bool useLine = !( line.substr(0,1) == "#" || line.substr(0,1) == "!");
    if ( useLine ) {
      instream.clear();          // Reset from possible previous errors
      instream.str(line);        // Use line as source of input
      lineNum++;
      if ( lineNum == 1 ) {
        instream >> lumiTarget;
      }
      else if ( lineNum == 2 ) {
        instream >> ns >> lumiSig;
	s = ns*lumiTarget/lumiSig;
      }
      else {
        double m_i, lumiMC_i, w_i;
        instream >> m_i >> lumiMC_i >> w_i;
        double tau_i = lumiMC_i / lumiTarget;
        double b_i = m_i / tau_i;
        b.push_back(b_i);
        m.push_back(m_i);
        tau.push_back(tau_i);
        w.push_back(w_i);
      }
    }
  }
  int numBck = m.size();

// Set up data n for a given mu, e.g., n = mu*s + btot

  double muData = 1.0;       // for creating data set
  double n = muData*s;
  for (int i=0; i<numBck; i++){
    //    cout << "i, w, b, m = " << i << "  " << w[i] << "  " << b[i] << "  " << m[i] << endl;
    n += w[i]*b[i];          // add best estimate of btot = sum_i w_i*b_i
  }
  // cout << "n = " << n << endl;

  double muTest = 0.0;       // for defining lambda(mu) (0 for discovery)
  double Z = getSignificance(muTest, n, s, m, tau, w);
  cout << "Z   = " << Z << endl;

  return 0;

}