void FeldmanCousins()
{
 // Example macro of using the TFeldmanCousins class in root.
 //
 // get a FeldmanCousins calculation object with the default limits
 // of calculating a 90% CL with the minimum signal value scanned 
 // = 0.0 and the maximum signal value scanned of 50.0
 //Author : Adrian John Bevan <bevan@SLAC.Stanford.EDU>
  
 if (!gROOT->GetClass("TFeldmanCousins")) gSystem->Load("libPhysics");
 
 TFeldmanCousins f;
 TGraph* FCdata       = new TGraph(20);
 TGraph* FCmc         = new TGraph(20);
 TGraph* FCdata_sys   = new TGraph(20);
 TGraph* FCmc_sys     = new TGraph(20);
 TGraph* FCdata_l     = new TGraph(20);
 TGraph* FCmc_l       = new TGraph(20);
 TGraph* FCdata_sys_l = new TGraph(20);
 TGraph* FCmc_sys_l   = new TGraph(20);
 
 // calculate either the upper or lower limit for 10 observerd
 // events with an estimated background of 3.  The calculation of
 // either upper or lower limit will return that limit and fill
 // data members with both the upper and lower limit for you.
 // Double_t Nobserved   = 129870.0;
 // Double_t Nbackground = 129870.0 + 303.0/1.994925;
 Double_t Nobserved   = 129870.0;
 Double_t Nbackground = 129870.0 + 303.0/(1.994925);
 f.SetMuStep(1);
 f.SetMuMax(155000);
 f.SetMuMin(0);

 
 cout << "For " <<  Nobserved << " data observed with and estimated background"<<endl;
 cout << "of " << Nbackground << " candidates, the Feldman-Cousins method of "<<endl;

 Double_t ul = f.CalculateUpperLimit(Nobserved, Nbackground);
 Double_t ll = f.GetLowerLimit();
 
 TMarker* m_data_stat = new TMarker(-303.0, ul * 1.994925, 20);

 cout << "calculating confidence limits gives:"<<endl;
 cout << "\tUpper Limit = " <<  ul*(1.994925) << endl;
 cout << "\tLower Limit = " <<  ll*(1.994925) << endl;
 cout << "at the 90% CL"<< endl;

 Nobserved   = 127772.;
 Nbackground = 127772. + 77./ 1.995938;
 // Double_t Nobserved   = 10;
 // Double_t Nbackground = 5;

 
 cout << "For " <<  Nobserved << " data observed with and estimated background"<<endl;
 cout << "of " << Nbackground << " candidates, the Feldman-Cousins method of "<<endl;

 ul = f.CalculateUpperLimit(Nobserved, Nbackground);
 ll = f.GetLowerLimit();

 TMarker* m_mc_stat = new TMarker(-77., ul * 1.995938, 21);

 cout << "calculating confidence limits gives:"<<endl;
 cout << "\tUpper Limit = " <<  ul*(1.995938) << endl;
 cout << "\tLower Limit = " <<  ll*(1.995938) << endl;
 cout << "at the 90% CL"<< endl;
 
 Nobserved   = 129870.0;
 // Nbackground = 129870.0*(1. - 0.0023*0.0023) + 303.0/1.994925;
 Nbackground = 129870.0 + 303.0/(2.6819);
 f.SetMuStep(1);
 f.SetMuMax(155000);
 f.SetMuMin(0);


 ul = f.CalculateUpperLimit(Nobserved, Nbackground);
 ll = f.GetLowerLimit();

 TMarker* m_data_sys = new TMarker(-303.0, ul * 2.6819, 22);
 
 cout << "For " <<  Nobserved << " data observed with and estimated background"<<endl;
 cout << "of " << Nbackground << " candidates, the Feldman-Cousins method of "<<endl;
 cout << "calculating confidence limits gives:"<<endl;
 cout << "\tUpper Limit = " <<  ul*(2.6819) << endl;
 cout << "\tLower Limit = " <<  ll*(2.6819) << endl;
 cout << "at the 90% CL"<< endl;

 Nobserved   = 127772;
 // Nbackground = 127772*(1. - 0.0037*0.0037) + 77/1.995938;
 Nbackground = 127772. + 77./3.7451;
 // Double_t Nobserved   = 10;
 // Double_t Nbackground = 5;

 ul = f.CalculateUpperLimit(Nobserved, Nbackground);
 ll = f.GetLowerLimit();
 TMarker* m_mc_sys = new TMarker(-77., ul * 3.7451, 23);

 cout << "For " <<  Nobserved << " data observed with and estimated background"<<endl;
 cout << "of " << Nbackground << " candidates, the Feldman-Cousins method of "<<endl;
 cout << "calculating confidence limits gives:"<<endl;
 cout << "\tUpper Limit = " <<  ul*3.7451 << endl;
 cout << "\tLower Limit = " <<  ll*3.7451 << endl;
 cout << "at the 90% CL"<< endl;


 for(int i=0; i<20; i++){
   double m = -1900. + 200. * i;
   Double_t Nobserved   = 129870.0;
   Double_t Nbackground = 129870.0 - m/(1.994925);
   ul = f.CalculateUpperLimit(Nobserved, Nbackground);
   ll = f.GetLowerLimit();
   FCdata->SetPoint(i, m, ul * 1.994925);
   FCdata_l->SetPoint(i, m, ll * 1.994925);
   
   Double_t Nbackground = 129870.0 - m/(2.6819);
   ul = f.CalculateUpperLimit(Nobserved, Nbackground);
   ll = f.GetLowerLimit();
   FCdata_sys->SetPoint(i, m, ul * 2.6819);
   FCdata_sys_l->SetPoint(i, m, ll * 2.6819);

   Double_t Nobserved   = 127772.;
   Double_t Nbackground = 127772. - m/(1.995938);
   ul = f.CalculateUpperLimit(Nobserved, Nbackground);
   ll = f.GetLowerLimit();
   FCmc->SetPoint(i, m, ul * 1.995938);
   FCmc_l->SetPoint(i, m, ll * 1.995938);
   
   Double_t Nbackground = 127772. - m/(3.7451);
   ul = f.CalculateUpperLimit(Nobserved, Nbackground);
   ll = f.GetLowerLimit();
   FCmc_sys->SetPoint(i, m, ul * 3.7451);
   FCmc_sys_l->SetPoint(i, m, ll * 3.7451);

 }
 
 TFile* file = new TFile("FCoutput.root","recreate");

 FCdata->SetLineColor(1);
 FCdata->SetLineStyle(1);
 FCdata_sys->SetLineColor(2);
 FCdata_sys->SetLineStyle(2);
 FCdata_l->SetLineColor(1);
 FCdata_l->SetLineStyle(1);
 FCdata_sys_l->SetLineColor(2);
 FCdata_sys_l->SetLineStyle(2);

 FCmc->SetLineColor(4);
 FCmc->SetLineStyle(4);
 FCmc_sys->SetLineColor(6);
 FCmc_sys->SetLineStyle(6);
 FCmc_l->SetLineColor(4);
 FCmc_l->SetLineStyle(4);
 FCmc_sys_l->SetLineColor(6);
 FCmc_sys_l->SetLineStyle(6);
 
 TLegend* leg = new TLegend(0.15, 0.6, 0.4, 0.85);
 leg->AddEntry(FCdata,"90% C.L., statistics only","l");
 leg->AddEntry(FCdata_sys,"90% C.L., systematics added", "l");
 leg->AddEntry(m_data_stat,"Measured 90% Upper Limit", "p");
 leg->AddEntry(m_data_sys,"Measured 90% Upper Limit", "p");
 leg->SetFillColor(10);

 TCanvas* c1 = new TCanvas();
 FCdata_sys->Draw("al");
 FCdata_sys->SetTitle(";Measured Number of Pions; True Number of Pions");
 FCdata_sys->GetHistogram()->GetXaxis()->SetTitleSize(0.05);
 FCdata_sys->GetHistogram()->GetYaxis()->SetTitleSize(0.05);
 FCdata_sys->GetHistogram()->GetXaxis()->SetLabelSize(0.05);
 FCdata_sys->GetHistogram()->GetYaxis()->SetLabelSize(0.05);
 FCdata_sys->GetHistogram()->GetXaxis()->SetRangeUser(-1900,1900);
 FCdata_sys->GetHistogram()->GetYaxis()->SetRangeUser(0,2500);
 FCdata->Draw("lsame");
 FCdata_sys_l->Draw("lsame");
 FCdata_l->Draw("lsame");
 m_data_stat->Draw("psame");
 m_data_sys->Draw("psame");
 leg->Draw("same");
 c1->Print("data_feldmancousins_limits.eps");
 
 c1->Clear();
 leg->Clear();

 
 leg->AddEntry(FCmc,"90% C.L., Statistics Only","l");
 leg->AddEntry(FCmc_sys,"90% C.L., Systematics Added", "l");
 leg->AddEntry(m_mc_stat,"Measured 90% Upper Limit Stat. Only", "p");
 leg->AddEntry(m_mc_sys,"Measured 90% Upper Limit Sys. Added", "p");

 FCmc_sys->Draw("al");
 FCmc_sys->SetTitle(";Measured Number of Pions; True Number of Pions");
 FCmc_sys->GetHistogram()->GetXaxis()->SetTitleSize(0.05);
 FCmc_sys->GetHistogram()->GetYaxis()->SetTitleSize(0.05);
 FCmc_sys->GetHistogram()->GetXaxis()->SetLabelSize(0.05);
 FCmc_sys->GetHistogram()->GetYaxis()->SetLabelSize(0.05);
 FCmc_sys->GetHistogram()->GetXaxis()->SetRangeUser(-1900,1900);
 FCmc_sys->GetHistogram()->GetYaxis()->SetRangeUser(0,2500);
 FCmc->Draw("lsame");
 FCmc_sys_l->Draw("lsame");
 FCmc_l->Draw("lsame");
 m_mc_stat->Draw("psame");
 m_mc_sys->Draw("psame");
 leg->Draw("same");
 c1->Print("mc_feldmancousins_limits.eps");

 file->cd();
 FCdata->Write();
 FCdata_sys->Write();
 FCmc->Write();
 FCmc_sys->Write();
 FCdata_l->Write();
 FCdata_sys_l->Write();
 FCmc_l->Write();
 FCmc_sys_l->Write();
 m_data_stat->Write();
 m_data_sys->Write();
 m_mc_stat->Write();
 m_mc_sys->Write();
 
 // f->Write();
 file->Close();

}