// ----------------------------------------------------------------------------------- //
/*
  ROOT macro illustrating simple example of Neyman-Pearson (Likelihood ratio) test:

  For more information, see:
    Barlow: Chapter 8 (following the example in 8.1.5)
    http://en.wikipedia.org/wiki/Likelihood-ratio_test

  Author: Troels C. Petersen (NBI/CERN)
  Email:  Troels.Petersen@cern.ch
  Date:   20th of September 2010
*/
// ----------------------------------------------------------------------------------- //

#include <TFile.h>
#include <TH1F.h>
#include <TRandom3.h>


//---------------------------------------------------------------------------------- 
double sqr(double a) {
  return a*a;
}


//---------------------------------------------------------------------------------- 
void NeymanPearson() 
//---------------------------------------------------------------------------------- 
{
  // Setting of general plotting style:
  gStyle->SetCanvasColor(0);
  gStyle->SetFillColor(0);
  // Setting what to be shown in statistics box:
  gStyle->SetOptStat("e");
  gStyle->SetOptFit(1111);

  // Set parameters:
  int Nopal   = 1000;
  int Nquartz = 1000;

  double mu_opal    = 2.2;
  double sig_opal   = 0.2;
  double mu_quartz  = 2.6;
  double sig_quartz = 0.2;

  TRandom3 r;

  // Make histograms:
  TH1F *Hist_Dens = new TH1F( "Hist_Dens", "Density distribution", 100, 1.4, 3.4);


//----------------------------------------------------------------------------------
// Loop over process:
//----------------------------------------------------------------------------------

  for (int i=0; i < Nopal; i++) {
    // Get random number from opal density distributed:
    double x_opal = r.Gaus(mu_opal, sig_opal);
    Hist_Dens->Fill(x_opal);
  }

  for (int i=0; i < Nquartz; i++) {
    // Get random number from quartz density distributed:
    double x_quartz = r.Gaus(mu_quartz, sig_quartz);
    Hist_Dens->Fill(x_quartz);
  }


//---------------------------------------------------------------------------------- 
// Plot histograms on screen:
//---------------------------------------------------------------------------------- 

  // Making a new window (canvas):
  // --------------------------------------------------------------------
  TCanvas* c0 = new TCanvas("c0", "", 100, 20, 600, 400);

  // Setting line width and color and axis titles of histograms:
  Hist_Dens->GetXaxis()->SetTitle("Density (g/cc)");
  Hist_Dens->GetYaxis()->SetTitle("Frequency");
  Hist_Dens->SetLineWidth(2);
  Hist_Dens->SetLineColor(4);

  // Drawing histograms with errors in the same plot:
  Hist_Dens->Draw();

  c0->Update();
  // c0->SaveAs("Dist_Density.png");




//---------------------------------------------------------------------------------- 
// Write histograms to file:
//---------------------------------------------------------------------------------- 
  /*
  TFile *file = new TFile("results.root","RECREATE","Histograms");

  Hist_Dens->Write();
  
  file->Write();
  file->Close();
  */
}



//---------------------------------------------------------------------------------- 
/*

Yes, once again, read through the program, and see if you understand what is going on!!!

Questions:
----------
 1) Look at the initial distribution - is it easy to see, that this is in fact two
    Gaussian distributions? Try to fit this distribution, first with a single Gaussian
    and then combining two. Remember to give reasonable initial values.

 2) Given you knowledge of the opal and quartz density and the measuring width (thus
    making it a simple hypothesis), what value would you cut at in density to get
    maximum 5% of quartz in your opal sample? Check in the simulation, if your value
    actually attained this?

 3) What would you set the cut at, if you wanted a factor 20 rejection against quartz
    (i.e. a power of 95%)?


Advanced questions:
-------------------
 1) Imagine that you had some other independent way of determining if it was opal or
    quartz, which had the same separation (i.e. two Gaussians with the same values).
    Try to simulate this situation, and see if you can pick out the opal without
    using your prior knowledge of the parameters, but simply that you know that opal
    is lighter than quartz (and similarly for the other new variable/method).

*/
//----------------------------------------------------------------------------------