void plot_sigback2(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 npts1 200;
   #define npts2 1500;

   // use flag to select weight

   Double_t wflag=0;    // weight = 1
   // Double_t wflag=1;    // weight = Gaussian

   // define histograms

   Int_t nbins=100;

   Double_t b1=4.;    // unis of 1/x
   Double_t b2=4.;    // unis of 1/x

   TH1F *expon1 = new TH1F("expon1","Projection Sig",nbins,0,10/b2);
   expon1->SetTitle("");
   TH2F *sig = new TH2F("sig","Signal",nbins,0,10/b2,nbins,0,10/b2);
   sig->SetTitle("");
   TH1F *integ1 = new TH1F("integ1","Integral Sig",nbins,0,10/b2);
   integ1->SetTitle("");

   TH1F *expon2 = new TH1F("expon2","Projection Back",nbins,0,10/b2);
   expon2->SetTitle("");
   TH2F *back = new TH2F("back","Background",nbins,0,10/b2,nbins,0,10/b2);
   back->SetTitle("");
   TH1F *integ2 = new TH1F("integ2","Integral Back",nbins,0,10/b2);
   integ2->SetTitle("");

   TH1F *sigb = new TH1F("sigb","Signal to Background above cut",nbins,0,10/b2);
   sigb->SetTitle("");
   TH1F *sigsqrtb = new TH1F("sigsqrtb","Signal to sqrt(Background) above cut",nbins,0,10/b2);
   sigsqrtb->SetTitle("");

   TF2 *sigf = new TF2("sigf",sigf,0,10/b2,0,10/b2,3);
   sigf->SetParameters(b1,b2,wflag);

   TF2 *backf = new TF2("backf",backf,0,10/b2,0,10/b2,3);
   backf->SetParameters(b1,b2,wflag);

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

   Double_t tk;
   Double_t tphi;

   for (j=0;j<npts1;j++) {
      sigf->GetRandom2(tk,tphi);
      sig->Fill(tk,tphi);
      expon1->Fill(tk);
      }

   for (j=0;j<npts2;j++) {
      backf->GetRandom2(tk,tphi);
      back->Fill(tk,tphi);
      expon2->Fill(tk);
      }

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

   c3->Divide(2,2);
   c3->cd(1);
   c3_1->SetBorderMode(0);
   c3_1->SetFillColor(0);
   c3_1->SetLogz();

   sig->SetTitle("Signal");
   // sig->GetXaxis()->SetRangeUser(1.0,2.5);
   // sig->GetYaxis()->SetRangeUser(ymin,ymax);
   sig->GetXaxis()->SetTitleSize(0.05);
   sig->GetYaxis()->SetTitleSize(0.05);
   sig->GetYaxis()->SetTitleOffset(1.5);
   sig->GetYaxis()->SetTitle("-t_{#phi}");
   sig->GetXaxis()->SetTitle("-t_{K}");
   sig->GetXaxis()->SetNdivisions(505);
   sig->SetMarkerColor(2);
   sig->SetMarkerStyle(21);
   sig->Draw("colz");
   // sigvs2->Draw("colz");

   c3->cd(2);
   c3_2->SetBorderMode(0);
   c3_2->SetFillColor(0);
   c3_2->SetLogz();

   back->SetTitle("Background");
   // back->GetXaxis()->SetRangeUser(1.0,2.5);
   // back->GetYaxis()->SetRangeUser(ymin,ymax);
   back->GetXaxis()->SetTitleSize(0.05);
   back->GetYaxis()->SetTitleSize(0.05);
   back->GetYaxis()->SetTitleOffset(1.5);
   back->GetYaxis()->SetTitle("-t_{#phi}");
   back->GetXaxis()->SetTitle("-t_{K}");
   back->GetXaxis()->SetNdivisions(505);
   back->SetMarkerColor(2);
   back->SetMarkerStyle(21);
   back->Draw("colz");

   c3->cd(3);
   c3_3->SetBorderMode(0);
   c3_3->SetFillColor(0);
   expon1->SetTitle("");
   // expon1->GetXaxis()->SetRangeUser(1.0,2.5);
   // expon1->GetYaxis()->SetRangeUser(ymin,ymax);
   expon1->GetXaxis()->SetTitleSize(0.05);
   expon1->GetYaxis()->SetTitleSize(0.05);
   expon1->GetYaxis()->SetTitleOffset(1.5);
   expon1->GetYaxis()->SetTitle("events");
   expon1->GetXaxis()->SetTitle("-t_{k}");
   expon1->GetXaxis()->SetNdivisions(505);
   expon1->SetMarkerColor(2);
   expon1->SetMarkerStyle(21);
   expon1->Fit("expo","","",0.4,2.5);
   expon1->Draw();

   c3->cd(4);
   c3_4->SetBorderMode(0);
   c3_4->SetFillColor(0);
   expon2->SetTitle("");
   // expon2->GetXaxis()->SetRangeUser(1.0,2.5);
   // expon2->GetYaxis()->SetRangeUser(ymin,ymax);
   expon2->GetXaxis()->SetTitleSize(0.05);
   expon2->GetYaxis()->SetTitleSize(0.05);
   expon2->GetYaxis()->SetTitleOffset(1.5);
   expon2->GetYaxis()->SetTitle("events");
   expon2->GetXaxis()->SetTitle("-t_{k}");
   expon2->GetXaxis()->SetNdivisions(505);
   expon2->SetMarkerColor(2);
   expon2->SetMarkerStyle(21);
   expon2->Draw();
   expon2->Fit("expo","","",0.4,2.5);

   //
   TCanvas *c2 = new TCanvas("c2","c2 plot_sigback2",200,10,700,700);
   c2->SetBorderMode(0);
   c2->SetFillColor(0);
   // c2->SetGridx();
   // c2->SetGridy();
   // c2->SetLogz();
   c2->Divide(2,2);

   // get integrals

   for (j=0;j<nbins;j++) {
       Float_t sum1 = sig->ProjectionX()->Integral(0,j);
       Float_t xbin = sig->ProjectionX()->GetBinCenter(j+1);
       Float_t ybin1 = sig->ProjectionX()->GetBinContent(j+1);
       integ1->Fill(xbin,sum1);
       Float_t sum2 = back->ProjectionX()->Integral(0,j);
       Float_t xbin = back->ProjectionX()->GetBinCenter(j+1);
       Float_t ybin2 = back->ProjectionX()->GetBinContent(j+1);
       integ2->Fill(xbin,sum2);

       if (sum2 > 0) {
          sigb->Fill(xbin,sum1/sum2);
          sigsqrtb->Fill(xbin,sum1/sqrt(sum2));
          }
       }

   Float_t xmin=0;
   Float_t xmax=2;
   Float_t ymin=1;
   Float_t ymax=npts2*5;

   c2->cd(1);
   c2_1->SetBorderMode(0);
   c2_1->SetFillColor(0);
   c2_1->SetLogy();

   // integ1->Scale(1/(Double_t) npts1);
   // integ2->Scale(1/(Double_t) npts2);
   integ1->SetTitle("");
   // integ1->GetXaxis()->SetRangeUser(1.0,2.5);
   integ1->GetYaxis()->SetRangeUser(ymin,ymax);
   integ1->GetXaxis()->SetTitleSize(0.05);
   integ1->GetYaxis()->SetTitleSize(0.05);
   integ1->GetYaxis()->SetTitleOffset(1.5);
   integ1->GetYaxis()->SetTitle("Integral below -t_{K} cut");
   integ1->GetXaxis()->SetTitle("-t_{k} selection cut");
   integ1->GetXaxis()->SetNdivisions(505);
   integ1->SetMarkerColor(2);
   integ1->SetMarkerStyle(21);
   integ1->SetLineColor(2);
   integ1->Draw("");

   integ2->SetTitle("");
   /* integ2->GetXaxis()->SetRangeUser(1.0,2.5);
   integ2->GetYaxis()->SetRangeUser(ymin,ymax);
   integ2->GetXaxis()->SetTitleSize(0.05);
   integ2->GetYaxis()->SetTitleSize(0.05);
   integ2->GetYaxis()->SetTitleOffset(1.5);
   integ2->GetYaxis()->SetTitle("Integral below -t_{K} cut");
   integ2->GetXaxis()->SetTitle("-t_{k} selection cut");
   integ2->GetXaxis()->SetNdivisions(505);
   integ2->SetMarkerColor(2);
   integ2->SetMarkerStyle(21);*/
   integ2->SetLineColor(2);
   integ2->SetLineColor(4);
   integ2->Draw("same");
   sprintf (string,"Signal, b1=%.1f GeV^{-1}\n",b1);
   t1 = new TLatex(0.25,0.7,string);
   t1->SetTextColor(2);
   t1->SetNDC();
   t1->Draw();
   sprintf (string,"Background, b2=%.1f GeV^{-1}\n",b2);
   t1 = new TLatex(0.25,0.65,string);
   t1->SetTextColor(4);
   t1->SetNDC();
   t1->Draw();
   if (wflag) {
      sprintf (string,"Weight=Gaussian");
      }
   else {
      sprintf (string,"Weight=none");
   }
   t1 = new TLatex(0.25,0.6,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();
   sprintf (string,"Signal events=%d\n",npts1);
   t1 = new TLatex(0.25,0.25,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();
   sprintf (string,"Background events=%d\n",npts2);
   t1 = new TLatex(0.25,0.2,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();

   Float_t ymin=1e-2;
   Float_t ymax=1;
   c2->cd(2);
   c2_2->SetBorderMode(0);
   c2_2->SetFillColor(0);
   c2_2->SetLogy();
   // sigb->GetXaxis()->SetRangeUser(xmin,xmax);
   sigb->GetYaxis()->SetRangeUser(ymin,ymax);
   sigb->SetTitle("");
   sigb->GetXaxis()->SetTitleSize(0.05);
   sigb->GetYaxis()->SetTitleSize(0.05);
   sigb->GetYaxis()->SetTitleOffset(1.5);
   sigb->GetYaxis()->SetTitle("S/B below -t_{K} cut");
   sigb->GetXaxis()->SetTitle("-t_{k} selection cut");
   sigb->GetXaxis()->SetNdivisions(505);
   sigb->SetMarkerColor(2);
   sigb->SetMarkerStyle(21);
   sigb->SetLineColor(2);
   sigb->Draw();
   sprintf (string,"Signal/Background\n");
   t1 = new TLatex(0.25,0.8,string);
   t1->SetTextColor(4);
   t1->SetNDC();
   t1->Draw();

   Float_t ymin=0;
   Float_t ymax=10;
   if (!wflag) {
       Float_t ymax=20;
       Float_t ymin=0;
   }
   c2->cd(3);
   c2_3->SetBorderMode(0);
   c2_3->SetFillColor(0);
   // c2_3->SetLogy();
   // sigsqrtb->GetXaxis()->SetRangeUser(xmin,xmax);
   sigsqrtb->GetYaxis()->SetRangeUser(ymin,ymax);
   sigsqrtb->SetTitle("");
   sigsqrtb->GetXaxis()->SetTitleSize(0.05);
   sigsqrtb->GetYaxis()->SetTitleSize(0.05);
   sigsqrtb->GetYaxis()->SetTitleOffset(1.5);
   sigsqrtb->GetYaxis()->SetTitle("S/#sqrt{B} below -t_{K} cut");
   sigsqrtb->GetXaxis()->SetTitle("-t_{k} selection cut");
   sigsqrtb->GetXaxis()->SetNdivisions(505);
   sigsqrtb->SetMarkerColor(2);
   sigsqrtb->SetMarkerStyle(21);
   sigsqrtb->SetLineColor(2);
   sigsqrtb->Draw();
   sprintf (string,"Signal/Sqrt(Background)\n");
   t1 = new TLatex(0.2,0.8,string);
   t1->SetTextColor(4);
   t1->SetNDC();
   t1->Draw();

   Int_t b1int = b1*10+0.5;
   Int_t b2int = b2*10+0.5;
   Int_t wflagint = wflag;

   sprintf(string,"plot_sigback2_c2_b1%d_b2%d_w%d.eps",b1int,b2int,wflagint);
   c2->SaveAs(string);
   sprintf(string,"plot_sigback2_c2_b1%d_b2%d_w%d.gif",b1int,b2int,wflagint);
   c2->SaveAs(string);

   sprintf(string,"plot_sigback2_c3_b1%d_b2%d_w%d.eps",b1int,b2int,wflagint);
   c3->SaveAs(string);
   sprintf(string,"plot_sigback2_c3_b1%d_b2%d_w%d.gif",b1int,b2int,wflagint);
   c3->SaveAs(string);

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

   Double_t b1=par[0];
   Double_t b2=par[1];
   Double_t wflag=par[2];
   Double_t x1=x[0];
   Double_t x2=x[1];

   // take y1 to be the Theta, y2 to be the phi

   char string[256];
   Double_t pi=3.14159;
   Double_t sigma=0.1;
   Double_t tphicut=0.15;

   sprintf(string," bslope=%f, %f\n",b1,b2);
   // printf("sigf: %s",string);

   Double_t y1 = exp(-b1*x1/2);
   Double_t y2 = exp(-b2*x2/2);

   // compute simple correllation between tphi and tks
   // from Apr 23 plot from Moskov of tphi vs tks, simple parameterization gives
   //      tphi = -0.375 + 1.5 tks -> tks = 0.25 + 0.67tphi
   //      sigma of distribution ~ 0.1 GeV-1

   Double_t mean = -0.375 + 1.5*x1;
   Double_t weight = exp(-(x2-mean)*(x2-mean)/(2*sigma*sigma))/(sqrt(2*pi)*sigma);
   if (x2 < tphicut) weight = 0;
   if (! wflag) Double_t weight=1;


   return weight*y1*y2;
}
Double_t backf (Double_t *x, Double_t *par) {

   Double_t b1=par[0];
   Double_t b2=par[1];
   Double_t wflag=par[2];
   Double_t x1=x[0];
   Double_t x2=x[1];

   // take y1 to be the Theta, y2 to be the phi

   char string[256];;
   Double_t pi=3.14159;
   Double_t sigma=0.1;
   Double_t tphicut=0.15;

   sprintf(string," bslope=%f, %f\n",b1,b2);
   // printf("backf: %s",string);

   Double_t y1 = exp(-b1*x1/2);
   Double_t y2 = exp(-b2*x2/2);

   // compute simple correllation between tphi and tks
   // from Apr 23 plot from Moskov of tphi vs tks, simple parameterization gives
   //      tphi = -0.375 + 1.5 tks -> tks = 0.25 + 0.67tphi
   //      sigma of distribution ~ 0.1 GeV-1

   Double_t mean = -0.375 + 1.5*x1;
   Double_t weight = exp(-(x2-mean)*(x2-mean)/(2*sigma*sigma))/(sqrt(2*pi)*sigma);
   if (x2 < tphicut) weight = 0;
   if (! wflag) Double_t weight=1;

   return weight*y2*y2;
}