//Test program for random number generators (spped and quality)
//The program get n random number pairs x and y i [0,1]
//It counts the ratio of pairs in the circle of diameter 1
//compared to the total number of pairs.
//This ratio must be Pi/4
//The test shows the graph of the difference of this ratio compared to PI.
//To test only the speed, call pirndmthread(1) (default)
//To test quality and speed call pirndmthread(50)
//    root pirndmthread.C+
//or
//    root "pirndmthread.C+(10)"
//
//Author: Rene Brun
      
#include "TROOT.h"
#include "TStopwatch.h"
#include "TMath.h"
#include "TRandom1.h"
#include "TRandom2.h"
#include "TRandom3.h"
#include "TCanvas.h"
#include "TH2.h"
#include "TGraph.h"
#include "TSystem.h"
#include "TLegend.h"
#include "TPaveLabel.h"
#include "TThread.h"
   
TLegend *legend = 0;
TCanvas *c1 = 0;
Double_t cputot = 0;

TThread *th1, *th2, *th3, *th4, *th5;

typedef struct {
   char name[64];
   TRandom *r;
   Long64_t n;
   Int_t color;
} args_t;

args_t args1, args2, args3, args4, args5;

//____________________________________________________________________
void *piRandom(void *ptr) 
{
   TStopwatch timer;
   args_t *args = (args_t *)ptr;
   timer.Start();
   TThread::Lock();
   TGraph *gr = new TGraph();
   TThread::UnLock();
   gr->SetMarkerStyle(20);
   gr->SetMarkerSize(0.7);
   gr->SetMarkerColor(args->color);
   gr->SetLineColor(args->color);
   gr->SetLineWidth(2);
   
   Int_t k = 0;   
   Double_t diffpi;
   Long64_t npi   = 0;
   Double_t pi = TMath::Pi();
   const Int_t NR = 20000;
   const Int_t NR2 = NR/2;
   Double_t rn[NR];
   Long64_t i = 0;
   while (i<=args->n) {
      i += NR2;
      args->r->RndmArray(NR,rn);
      for (Int_t j=0;j<NR;j+=2) {
         if (rn[j]*rn[j]+rn[j+1]*rn[j+1] <= 1) npi++;
      }
      if (i && i %100000000 == 0) {
         Double_t norm = 4./Double_t(i);
         diffpi   = norm*npi - pi;
         gr->SetPoint(k,i,diffpi);
         if (k ==0) {
            TThread::Lock();
            gr->Draw("lp");
            TThread::UnLock();
         }
         c1->Modified();
         c1->Update();
         k++;
      }
      // yield execution to another thread that is ready to run
      // MSDN :
      // Sleep(0) causes the thread to relinquish the remainder of its 
      // time slice to any other thread of equal priority that is ready to run. 
      // If there are no other threads of equal priority ready to run, 
      // the function returns immediately, and the thread continues execution.
      gSystem->Sleep(0);
   }
   timer.Stop();
   Double_t cpu = timer.CpuTime();
   cputot += cpu;
   Double_t nanos = 1.e9*cpu/Double_t(args->n);
   legend->AddEntry(gr,Form("%-14s: %6.1f ns/call",args->name,nanos),"lp");
   c1->Modified();
   c1->Update();
   TThread::Printf("RANDOM = %s : RT=%7.3f s, Cpu=%7.3f s\n",args->name,timer.RealTime(),cpu);
   return 0;
}
         
//________________________________________________________________________
void ErrorBand(Long64_t n) 
{
   Int_t np = 40;
   TGraph *g = new TGraph(2*np+2);
   Double_t xmax = Double_t(n)/Double_t(np);
   for (Int_t i=1;i<=np;i++) {
      Double_t x = i*xmax;
      Double_t e = 1./TMath::Sqrt(x);
      g->SetPoint(i,x,e);
      g->SetPoint(2*np-i+1,x,-e);
   }
   Double_t x0 = 0.1*xmax;
   Double_t e0 = 1./TMath::Sqrt(x0);
   g->SetPoint(0,x0,e0);
   g->SetPoint(2*np+1,0,-e0);
   g->SetPoint(2*np+2,0,e0);
   g->SetFillColor(1);
   g->SetFillStyle(3002);
   g->Draw("f");
}

// thread to run Ps(): perform every 5 seconds a TThread::Ps()
void *top(void *)
{
   TThread::Printf("Start of top");

   for (Int_t i = 0; ;i++) {
      if ((th1->GetState() != TThread::kRunningState) &&
          (th2->GetState() != TThread::kRunningState) &&
          (th3->GetState() != TThread::kRunningState) &&
          (th4->GetState() != TThread::kRunningState) &&
          (th5->GetState() != TThread::kRunningState))
         break;
      gSystem->Sleep(100);
   }
   TThread::Printf("End of top");

   Double_t reftime = 9.25 + 12.01 + 11.91 + 31.0 + 57.02;
   const Double_t rootmarks = 900*Double_t(5.0)*reftime/cputot;
   TPaveLabel *pl = new TPaveLabel(0.2,0.92,0.78,0.98,Form("cpu time = %6.1fs - rootmarks = %6.1f",cputot,rootmarks),"brNDC");
   pl->Draw();  
   printf("******************************************************************\n");
   printf("*  ROOTMARKS =%6.1f   *  Root%-8s  %d/%d\n",rootmarks,gROOT->GetVersion(),
         gROOT->GetVersionDate(),gROOT->GetVersionTime());
   printf("******************************************************************\n");

   printf("Time at the end of job = %f seconds\n",cputot);
   c1->Print("pirndmthread.root");
   c1->Print("pirndmthread.gif");
   return 0;
}

//________________________________________________________________________
void pirndmthread(Long64_t n1=1) 
{
   Long64_t n = n1*200000000;
   c1 = new TCanvas("c1");
   c1->SetLeftMargin(0.12);
   c1->SetFrameFillColor(41);
   c1->SetFrameBorderSize(6);
   c1->SetGrid();
   Double_t dy = 1.5e-4;
   if (n1 < 4) dy *=2;
   TH2F *frame = new TH2F("h","",100,0,1.1*n,100,-dy,dy);
   frame->GetXaxis()->SetTitle("Number of TRandom calls");
   frame->GetYaxis()->SetTitle("Difference with #pi");
   frame->GetYaxis()->SetTitleOffset(1.3);
   frame->GetYaxis()->SetDecimals();
   frame->SetStats(0);
   frame->Draw();
   legend = new TLegend(0.6,0.7,0.88,0.88);
   legend->Draw();
   
   ErrorBand(n);

   sprintf(args1.name,"TRandom()");
   args1.r = new TRandom();
   args1.n = n;
   args1.color = kRed;
   printf("Starting Thread 1\n");
   th1 = new TThread("th1", piRandom, (void*) &args1);
   th1->Run();

   sprintf(args2.name,"TRandom1()");
   args2.r = new TRandom1();
   args2.n = n;
   args2.color = kBlack;
   printf("Starting Thread 2\n");
   th2 = new TThread("th2", piRandom, (void*) &args2);
   th2->Run();

   sprintf(args3.name,"TRandom2()");
   args3.r = new TRandom2();
   args3.n = n;
   args3.color = kMagenta;
   printf("Starting Thread 3\n");
   th3 = new TThread("th3", piRandom, (void*) &args3);
   th3->Run();

   sprintf(args4.name,"TRandom3()");
   args4.r = new TRandom3();
   args4.n = n;
   args4.color = kBlue;
   printf("Starting Thread 4\n");
   th4 = new TThread("th4", piRandom, (void*) &args4);
   th4->Run();

   sprintf(args5.name,"TRandom3(0)");
   args5.r = new TRandom3(0);
   args5.n = n;
   args5.color = kGreen;
   printf("Starting Thread 5\n");
   th5 = new TThread("th5", piRandom, (void*) &args5);
   th5->Run();

   // top thread
   TThread *thp = new TThread("top",top);
   thp->Run();

   TThread::Ps();
}