// ----------------------------------------------------------------------------------- //
/*
  ROOT macro for estimating the value of pi using a Monte Carlo technique.

  The program picks random points in the area [-1,1] x [-1,1], and determines which
  fraction of these are within the unit circle. This in turn gives a measure of pi
  with associated statistical uncertainty. Performing such "experiments" many times
  not only gives the value of pi, but also a feel for the use of Monte Carlo, and
  experience in calculating averages, RMSs, and the error on these.

  For more information see:
    P. R. Bevington: page 75-78

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

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



//---------------------------------------------------------------------------------- 
void PiEstimate() 
//---------------------------------------------------------------------------------- 
{
  TRandom3 r;
 
 // Set parameters:
  int Nexperiment = 1;          // Number of "experiments" determining pi
  int Npoints     = 1000;       // Number of points in determining pi

  double f, ef;
  double pi_estm, pi_error;
  double pi_true = 3.14159265358979323846264338328;

  double sum[3] = {0.0, 0.0, 0.0};

  // Make histograms:
  TH1F *Hist_PiDist  = new TH1F("Hist_PiDist",  "Hist_PiDist",   50, 2.9, 3.4);
  TH2F *Hist_HitDist = new TH2F("Hist_HitDist", "Hist_HitDist", 200, -1.0, 1.0, 200, -1.0, 1.0);


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

  for (int iexp=0; iexp < Nexperiment; iexp++) {

    int Nhit = 0;
    for (int j=0; j < Npoints; j++) {

      // Get two random numbers uniformly distributed in [-1,1]:
      double x1 = 2.0 * (r.Uniform() - 0.5);
      double x2 = 2.0 * (r.Uniform() - 0.5);

      // Plot 2D distribution of (x1,x2) for the first experiment:
      if (iexp == 0) Hist_HitDist->Fill(x1,x2);

      // Test if (x1,x2) is within unit circle:
      if (x1*x1 + x2*x2 < 1.0) Nhit++;
    }

    // Calculate the fraction of points within the circle and its error:
    f  = double(Nhit) / double(Npoints);
    ef = 0.0;   // This you have to calculate yourself!!! The formula is worth remembering!

    // From this we can get pi and its error:
    pi_estm  = 4.0 * f;
    pi_error = 4.0 * ef;

    // Plot these values and calculate mean and RMS:
    Hist_PiDist->Fill(pi_estm);
    sum[0] += 1.0;
    sum[1] += pi_estm;
    sum[2] += pi_estm * pi_estm;

    // Print first couple of pi measurements:
    if (iexp < 5)
      printf("  %2d. pi estimate:   %7.4f +- %6.4f \n", iexp+1, pi_estm, pi_error);
  }

  // Calculate average and RMS, and print the result compared to the true value of pi...

  // ...yourself!!!




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

  Hist_PiDist->Write();
  Hist_HitDist->Write();

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



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

First acquaint yourself with the program, and make sure that you understand what the
parameters "Nexperiment" and "Npoints" refere to! Also, before running the program,
calculate what precision you expect on pi, when using the number of points chosen in
the program.

Then, run the program, and then take a look at the result. Did you estimate match
the number given by the program?


Questions:
----------
 1) Try to run 100 experiments with 1000 points in each. What is the approximate
    uncertainty on pi in each experiment? What is the uncertainty on the AVERAGE
    of all 100 experiments?

 2) How do you expect the values of pi to distribute themselves? And is this the
    case in reality?

 3) Does it make any difference on the precision of the final pi value, whether
    you make many experiments with few points, or one experiment with many points,
    as long as the product of Nexperiment x Npoints remains constant?
    What could be the reason for choosing the former approach?

 4) Now try to use this method in three dimensions to estimate the constant in
    front of the pi*r^3 expression for the volume. Do you get 3/4?
    Increase the dimensionality (say up to 10), and see if you can figure out
    the constants needed to calculate the hyper-volumes!


Advanced questions:
-------------------
 1) How quickly is the convergence of the pi precision, and how does this compare
    to an algorithm, which does not use random numbers, but simply does a numerical
    integral over the region?


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