// ntotal = Number of events per trial

void SimALL(int ntotal = 300, double& sigma, double& error)
{
  const double all_in = 0.2;	// Input A_LL
  const double pb1 = 0.7;	// Polarization of beam 1
  const double pb2 = 0.7;	// Polarization of beam 2
  const double neff = 0.01;	// Average number of interactions for each unpolarized bunch crossing
  const double mu1 = (1+pb1*pb2*all_in)*neff; // Average yield of events per bunch crossing with polarization state "+"
  const double mu2 = (1-pb1*pb2*all_in)*neff; // Average yield of events per bunch crossing with polarization state "-"
  const int ntrial = 500;	// Number of independent simulations of asymmetry

  // Create A_LL histogram
  TH1* hAll = new TH1F("hAll",";#frac{1}{P_{b_{1}}P_{b_{2}}}#frac{N_{+}-N_{-}}{N_{+}+N_{-}};Number of trials",50,-1,1);

  // Run  ntrial independent simulations each with ntotal events
  for (int i = 0; i < ntrial; ++i) {
    int n1 = 0;			// Number of events with beam helicity "+"
    int n2 = 0;			// Number of events with beam helicity "-"
    while (n1 + n2 < ntotal) {
      if (gRandom->Poisson(mu1) > 0) ++n1;
      if (gRandom->Poisson(mu2) > 0) ++n2;
    }

    // Calculate A_LL
    double all = 1/(pb1*pb2)*(n1-n2)/(n1+n2);

    // Fill histogram
    hAll->Fill(all);
  }

  // Plot
  hAll->GetXaxis()->CenterTitle(true);
  hAll->GetYaxis()->CenterTitle(true);
  hAll->SetTitleOffset(1.5);
  hAll->Fit("gaus","W");

  TF1* gaus = hAll->GetFunction("gaus");
  gaus->SetLineColor(kRed);

  sigma = gaus->GetParameter(2);
  error = gaus->GetParError(2);

  delete hAll;
}

void SimALL2()
{
  gROOT->SetStyle("Plain");

  TCanvas* c1 = new TCanvas("c1","c1",500,500);
  c1->SetLeftMargin(0.2);

  gStyle->SetOptStat(10);
  gStyle->SetOptFit(111);

  TH1* hSigma = new TH1F("hSigma",";Number of events per trial;Gaussian #sigma of #frac{1}{P_{b_{1}}P_{b_{2}}}#frac{N_{+}-N_{-}}{N_{+}+N_{-}} distribution",700,0,700);

  for (int ntrial = 50; ntrial < 700; ntrial += 50) {
    double sigma;
    double error;
    SimALL(ntrial, sigma, error);
    hSigma->SetBinContent(ntrial+1,sigma);
    hSigma->SetBinError(ntrial+1,error);
  }

  hSigma->SetStats(false);
  hSigma->GetXaxis()->CenterTitle(true);
  hSigma->GetYaxis()->CenterTitle(true);
  hSigma->GetYaxis()->SetTitleOffset(2.0);
  hSigma->SetMinimum(0);
  hSigma->SetMaximum(0.35);
  hSigma->SetMarkerStyle(kOpenCircle);
  hSigma->Draw();

  TF1* f1 = new TF1("f1","1/(0.7*0.7)*1/sqrt(x)",0,700);
  f1->SetLineColor(kRed);
  f1->Draw("same");

  c1->Print("Sigma.png");
  c1->Print("Sigma.eps");
}