float buffon(int N, float needleLength){
  // This should implement the buffon method to estimate Pi.
  //---------------The following is deleted in the template ---------------------------------------------------
  float counter = 0;
  //N will be a large number, so book variables in advance
  float needlePos   = 0.;
  float needleAngle = 0.;
  float distX       = 0.;
  for (int ii = 0; ii < N; ii++){
    needlePos   = gRandom->Uniform(0,1);
    needleAngle = gRandom->Uniform(0,TMath::Pi()/2.);

    distX = TMath::Min(needlePos, 1-needlePos);

    if(TMath::Cos(needleAngle) * needleLength/2. > distX){
      counter += 1;
    }
  }
  return float(N/counter * 2. * needleLength);
  //return 3.14
}


float alternative(int N){
  // This should implement the alternative Version.
  // -------------------The following is deleted in the template -------------------------------------------
  float counter = 0.;
  for (int ii = 0; ii < N; ii++){
    if (pow(gRandom->Uniform(0,1),2) + pow(gRandom->Uniform(0,1),2) < 1){
      counter += 1.;
    }
  }
  return (counter/N) * 4;
  //return 3.14
}


int ex12 (){
  //First we show, that the principle is working 
  int N = 100000;
  float needleLength = 1.;

  cout << "Proof of Principle." << endl;
  cout << "Pi is : " << endl;
  cout << TMath::Pi() << endl;

  cout << "Buffon Method says: " << endl;
  cout << buffon (N, needleLength)<< endl;

  cout << "Alternative Method says: " << endl;
  cout << alternative(N) << endl;

  //Now lets make some double logarithmic plot to guess the speed of conversion.
  //Unfortunately it turns out the differences are quite small, except for extreme needle length and one method is clearly better than the other...
  static const int lengthSteps = 5;
  static const int ordersOfMag = 40;
  float xList[lengthSteps+1][ordersOfMag];
  float yList[lengthSteps+1][ordersOfMag];
  for (int ii = 0; ii < lengthSteps; ii++){
    for (int jj = 0; jj < ordersOfMag; jj++){
      xList[ii][jj] = int(10 * pow(1.2, (jj+1)));
      yList[ii][jj] = fabs(buffon(int(10 * pow(1.2, (jj+1))), 1. - 0.2* ii)- TMath::Pi());
    }
  }

  for (int jj = 0; jj < ordersOfMag; jj++){
    xList[lengthSteps][jj] = int(10 * pow(1.2, (jj+1)));
    yList[lengthSteps][jj] = fabs(alternative(int(10 * pow(1.2, (jj+1))))- TMath::Pi());
  }

  TCanvas mycan;
  mycan.cd();
  mycan.SetLogx();
  mycan.SetLogy();

  TGraph* myGraphArray [lengthSteps + 1];
  for (int ii = 0; ii < lengthSteps + 1; ii++){
    myGraphArray[ii] = new TGraph(ordersOfMag, xList[ii], yList[ii]);
    myGraphArray[ii]->SetLineColor(ii+1);
    myGraphArray[ii]->SetLineWidth(6);
    myGraphArray[ii]->SetMarkerStyle(ii);
    myGraphArray[ii]->SetMaximum(10);
    myGraphArray[ii]->SetMinimum(0.00001);
    if(ii ==0) {
      myGraphArray[ii]->Draw("AC*");
    }
    else {myGraphArray[ii]->Draw("same");
    }
  }
  mycan.SaveAs("Graph.root");

  return 0;

}