// ntotal = Number of events per trial

void SimALL(int ntotal = 300)
{
  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
  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(10);
  gStyle->SetOptFit(111);

  TCanvas* c1 = new TCanvas("c1","c1",500,500);
  c1->SetBottomMargin(0.15);

  hAll->GetXaxis()->CenterTitle(true);
  hAll->GetYaxis()->CenterTitle(true);
  hAll->GetXaxis()->SetTitleOffset(1.5);
  hAll->Fit("gaus","W");
  hAll->GetFunction("gaus")->SetLineColor(kRed);

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