void test_random(void)
{
// 
// plot the number of fits to a cdc segment that result in prob > prob_cut (currently set to 0.01).
//

#include <TMath.h>
#include <TRandom.h>

gROOT->Reset();
//TTree *Bfield = (TTree *) gROOT->FindObject("Bfield");
gStyle->SetPalette(1,0);
gStyle->SetOptStat(kTRUE);
gStyle->SetOptFit(kTRUE);
// gStyle->SetOptFit(1111);
gStyle->SetPadRightMargin(0.15);
gStyle->SetPadLeftMargin(0.15);
gStyle->SetPadBottomMargin(0.15);
gStyle->SetFillColor(0);
//

   char string[256];
   char filename[80];
   Int_t j,jj;
   
   #define npts 100000;

   // define histograms

   nbins=100;
   TH1F *uniform = new TH1F("uniform","uniform distribution",nbins,0,1);
   TH1F *decay = new TH1F("decay","exponential distribution",nbins,0,40);
   TH1F *bell = new TH1F("bell","bell distribution",nbins,-10,10);
   TH1F *counts= new TH1F("counts","Poisson distribution",nbins/10,0,10);

   //
   TCanvas *c1 = new TCanvas("c1","c1 Test random",200,10,700,700);
   c1->SetBorderMode(0);
   c1->SetFillColor(0);
   // c1->SetGridx();
   // c1->SetGridy();
   // c1->SetLogy();

   Double_t xmin=-1;
   Double_t xmax=1;

   // plot pulse function

   Double_t tau=4.;    // unis of 1/x
   Double_t mean=3;    // unis of 1/x
   Double_t sigma=2;     // units of x 
   Double_t Pmean=3;   // Poisson mean

   TRandom1 *r = new TRandom1();

   for (j=0;j<npts;j++) {
       Double_t x = r->Uniform();
       Double_t y = r->Exp(tau);
       Double_t z = r->Gaus(mean,sigma);
       Int_t n = r->Gaus(Pmean);

       uniform->Fill(x);
       decay->Fill(y);
       bell->Fill(z);
       counts->Fill((Double_t) n);
       }

   c1->Divide(2,2);
   c1->cd(1);
   c1_1->SetBorderMode(0);
   c1_1->SetFillColor(0);;
   uniform->GetYaxis()->SetRangeUser(0,2*npts/nbins);
   uniform->Fit("pol1");
   uniform->Draw();

   c1->cd(2);
   c1_2->SetBorderMode(0);
   c1_2->SetFillColor(0);
   decay->Fit("expo");
   decay->Draw();

   c1->cd(3);
   c1_3->SetBorderMode(0);
   c1_3->SetFillColor(0);
   bell->Fit("gaus");
   bell->Draw();

   c1->cd(4);
   c1_3->SetBorderMode(0);
   c1_3->SetFillColor(0);
   counts->Fit("gaus");
   counts->Draw();

   //
   TCanvas *c2 = new TCanvas("c2","c2 Test random",200,10,700,700);
   c2->SetBorderMode(0);
   c2->SetFillColor(0);
   // c2->SetGridx();
   // c2->SetGridy();
   // c2->SetLogy();

   // define input masses and kinematics

   Double_t mt = 0.938272;   // target units are GeV
   Double_t mb = 0;          // beam
   Double_t m1 = 0.938272;   // product 1
   Double_t m2 = 0.497614;   // product 2
   Double_t m3 = 0.497614;   // product 3
   Double_t Eb = 3;         // beam energy
   Double_t s = mt*mt + mb*mb + 2*Eb*mt; // total c.m. energy squared

   xmin = (m1+m2)*(m1+m2);
   xmax = (sqrt(s)-m3)*(sqrt(s)-m3);

   sprintf (string,"s=%f\nProduct 1 mass=%f\nProduct 2 mass=%f\nProduct 3 mass=%f\n",s,m1,m2,m3);
   printf ("test_random: %s",string);

   TF1 *Pm12 = new TF1("Pm12",Pm12,sqrt(xmin),sqrt(xmax),5);
   Pm12->SetParameters(s,m1,m2,m3,1);
   Double_t ymin = 0;
   Double_t ymax = 300;
   // Pm12->GetYaxis()->SetRangeUser(ymin,ymax);
   Pm12->DrawCopy();

   //
   TCanvas *c3 = new TCanvas("c3","c3 Test random",200,10,700,700);
   c3->SetBorderMode(0);
   c3->SetFillColor(0);
   //c3->SetGridx();
   //c3->SetGridy();
   //c3->SetLogy();

   TH1F *m12 = new TH1F("m12","m12",nbins,sqrt(xmin),sqrt(xmax));

   for (j=0;j<npts;j++) {
      Double_t rm12 = Pm12->GetRandom();
      m12->Fill(rm12);
      }
   Pm12->FixParameter(0,s);
   Pm12->FixParameter(1,m1);
   Pm12->FixParameter(2,m2);
   Pm12->FixParameter(3,m3);
   m12->Fit("Pm12","","",sqrt(xmin),sqrt(xmax));

   m12->SetTitle("");
   // m12->GetXaxis()->SetRangeUser(1.0,2.5);
   // m12->GetYaxis()->SetRangeUser(ymin,ymax);
   m12->GetXaxis()->SetTitleSize(0.05);
   m12->GetYaxis()->SetTitleSize(0.05);
   m12->GetYaxis()->SetTitleOffset(1.5);
   m12->GetYaxis()->SetTitle("events");
   m12->GetXaxis()->SetTitle("m_{pK}");
   m12->GetXaxis()->SetNdivisions(505);
   m12->SetMarkerColor(2);
   m12->SetMarkerStyle(21);
   m12->Draw();
      

}
Double_t Pm12 (Double_t *x, Double_t *par) {

   // generate phase space distribution in m12s
   // See Hagedorn Eq. 7-43

   char string[256];

   Double_t mstars=par[0];
   Double_t m1=par[1];
   Double_t m2=par[2];
   Double_t m3=par[3];
   Double_t norm=par[4];
   Double_t m12=x[0];
   Double_t mstar=sqrt(mstars);

   sprintf(string,"m1=%f, m2=%f, m3=%f, m12=%f, mstar=%f\n",m1,m2,m3,m12,mstar);
   // printf ("Pm12: %s",string);

   Double_t xin[1];
   Double_t parin[2];

   xin[0] = mstar;
   parin[0] = m12;
   parin[1] = m3;
   Double_t pm12 = sqrt(p2mstar(xin,parin));

   xin[0] = m12;
   parin[0] = m1;
   parin[1] = m2;
   Double_t pmj = sqrt(p2mstar(xin,parin));

   sprintf(string,"pm12=%f, pmj=%f\n",pm12,pmj);
   // printf ("Pm12: %s",string);

   // Double_t distr = pm12 * pmj / m12;  // distribution in m12s

   Double_t distr = pm12 * pmj; //  distribution in m12

   return norm*distr;

}
Double_t p2mstar (Double_t *x, Double_t *par) {

   // See Hagedorn Eq. 7-43

   Double_t m1=par[0];
   Double_t m2=par[1];
   Double_t mstar=x[0];

   Double_t mstars=mstar*mstar;
   Double_t m1s=m1*m1;
   Double_t m2s=m2*m2;

   char string[256];

   if (mstar < m1+m2) {
      // sprintf (string,"mstar=%f, m1=%f, m2=%f\n",mstar,m1,m2);
      // printf ("string=%s",string);
      return -1;
      }

   Double_t p2 = (mstars-(m1+m2)*(m1+m2))*(mstars - (m1-m2)*(m1-m2))/(4*mstars);
   
   sprintf (string,"p2mstar: m1=%f, m2=%f, mstar=%f, p2=%f\n",m1,m2,mstar, p2);
   // printf ("p2mstar: %s",string);

   return p2;
}