// File:  getDiscoverySignificance.cc
// Glen Cowan
// RHUL Physics

// Returns significance Z with which hypothesized value of mu=0 is rejected; 
// Z = Phi^-1(1 - p), p = p-value of mu, Phi^-1 = quantile of standard Gaussian. 
// 
// Inputs:     
// n     number of events seen in data, e.g. generate from Poisson(mu*s+b)
// s     expected number of signal events
// m     vector of numbers of events seen in subsidiary measurements.
// S1    sum of weights, w
// S2    sum of squares of w
// tau   defined by m_i ~ Poisson(tau_i*b/omega), where omega = E[w]

#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include <cstdlib>
#include <cmath>

#include "SigCalc.h"

double getDiscoverySignificance(double n, double s,
       vector<double> m, vector<double> tau, 
       vector<double> S1, vector<double> S2, int option){

  SigCalc* sc = new SigCalc (n, s, m, tau, S1, S2, option);

  // start values for fit now set from inside SigCalc

  /*  

  double sigma2Hat_i, lambdaHat_i, bHat_i, omegaHat_i;
  int numBkg = m.size();
  double btotHat = 0.;
  for (int i=0; i<numBkg; i++){ 

    double bHat_i = 0;
    if ( option < 100 ) {       // w ~ log-normal
      double lambdaHat_i;
      double sigma2Hat_i;
      if ( m[i] == 0 ) {
        sigma2Hat_i = 1.;
        lambdaHat_i = 1.;
      }
      else if ( m[i] == 1 ) {
        lambdaHat_i = S1[i]/m[i]; 
        double lambda0 = 0;              // nominal ln w
        sigma2Hat_i = (lambdaHat_i - lambda0)*(lambdaHat_i - lambda0);
      }
      else {
        lambdaHat_i = S1[i]/m[i]; 
        sigma2Hat_i = 
            (1./(m[i]-1)) * (S2[i] - S1[i]*S1[i]/m[i]);
      }
      double omegaHat_i = exp(lambdaHat_i + 0.5*sigma2Hat_i);
      bHat_i = m[i]*omegaHat_i/tau[i];
    }
    else if ( option >= 100 && option < 200 ) {     // w ~ Gaussian
      bHat_i = S1[i]/tau[i];
    }
    else if ( option >= 200 && option < 300 ) {     // bHat ~ Gaussian
      bHat_i = S1[i]/tau[i];
    }
    btotHat += bHat_i;
    cout << "prelim. i,  bHat_i = " << i << "  " << bHat_i << endl;
  }

  double muStart = (n - btotHat) / s;

 */

  double muHat;
  vector<double> bHat;
  vector<double> bHatHat;
  vector<double> lambdaHat;
  vector<double> lambdaHatHat;
  vector<double> sigmaHat;
  vector<double> sigmaHatHat;
  double q0Val = sc->q0(muHat, bHat, bHatHat, lambdaHat, 
			lambdaHatHat, sigmaHat, sigmaHatHat);
  double Z = sqrt(q0Val);


  // KLUDGE -- debug section

  //  cout << "n = " << sc->n() << endl;
  //  cout << "m = " << sc->m(0) << endl;
  //  cout << "S1 = " << sc->S1(0) << endl;
  //  cout << "S2 = " << sc->S2(0) << endl;

  /*
  while ( true ){

    double mu, b, lambda, sigma;
    vector<double> bVec(1);
    vector<double> lambdaVec(1);
    vector<double> sigmaVec(1);
    
    cout << "enter mu: ";
    cin >> mu;
    cout << "enter b: ";
    cin >> b;
    cout << "enter lambda "; 
    cin >> lambda;
    cout << "enter sigma ";
    cin >> sigma;
    bVec[0] = b;
    lambdaVec[0] = lambda;
    sigmaVec[0] = sigma;
    double logL = sc->lnL(mu, bVec, lambdaVec, sigmaVec); 
    cout << "logL = " << logL << endl;

  }

  */

  delete sc;
  return Z;

}