// 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 expected number of signal
// events s and the integrated luminosity L_data for which one wants 
// the significance (note this is in general not the effective luminosity
// of MC sample used to determine s; that quantity does not enter here).
//
// On subsequent lines the user lists for each background component the
// expected number of background events b corresponding to the luminosity 
// L_data and the effective luminosity of the MC sample L_MC used to generate 
// the background sample.  That is, the number of background events 
// accepted in the MC sample is m = b*L_MC/L_data.  The units for L_data
// and L_MC must be the same but are otherwise arbitrary.

#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 s, lumiData;
  vector<double> m;
  vector<double> b;
  vector<double> tau;
  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 >> s >> lumiData;
      }
      else {
        double b_i, lumiMC_i;
        instream >> b_i >> lumiMC_i;
        double tau_i = lumiMC_i / lumiData;
        double m_i = b_i * tau_i;
        b.push_back(b_i);
        m.push_back(m_i);
        tau.push_back(tau_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++){
    n += b[i];            // add best estimate of btot = sum_i b_i
  }

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

  return 0;

}