void gen_phithe(void)
{
// 
// Generate 1-dim background with gaussian signal superimposed to study significance 
//    as a function of various parameters.
//

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

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

   char string[256];
   char filename[80];
   char dataname[80];
   Int_t j,jj;
   
   Int_t nevents=6000;
   Double_t phi_norm=0.5;
   Double_t theta_norm=1.0;

   // define histograms

   Int_t nbins=100;
   Double_t xmin=1.4;
   Double_t xmax=2.0;
   Double_t binsize = (xmax-xmin)/nbins;

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

   // define input parameters

   Double_t mu=5.7;            // unis of 1/x
   Double_t sigma=0.035;      // units of x  
   Double_t gain=17;         // arbitrary
   Double_t xshift=1.47;     // x translation
   Double_t mean=1.54;       // gauss mean  
   Double_t sig=0.005;       // gauss sig  
   Double_t gnorm=40*sqrt(2*3.14159)*sig;       // gauss normalization

   sprintf (string,"s=%f\nmu=%f\nsigma=%f\ngain=%f\nxhift=%f\nmean=%f\nsig=%f\ngnorm=%f\n",mu,sigma,gain,xshift,mean,sig,gnorm);
   printf ("gen_phithe: %s",string);

   TF1 *PMMks = new TF1("PMMks",PMMks,xmin,xmax,7); 
   PMMks->SetParameters(mu,sigma,gain,xshift,gnorm,mean,sig);

   // get ratio of integral with nominal parameters to modified function.

   Double_t norm_integ = PMMks->Integral(xmin,xmax);
   gain=gain*phi_norm*phi_norm;         
   gnorm = gnorm*phi_norm*theta_norm;
   PMMks->SetParameters(mu,sigma,gain,xshift,gnorm,mean,sig);
   Double_t actual_integ = PMMks->Integral(xmin,xmax);
   Double_t ratio_integ = actual_integ/norm_integ;

   // scale number of generated events by ratio of integrals

   nevents = nevents * ratio_integ;
   printf ("phi_norm=%f, theta_norm=%f\n",phi_norm,theta_norm);
   printf ("Normalization integral=%f, Actual integral=%f, Ratio=%f, Nevents=%d\n",norm_integ,actual_integ,ratio_integ,nevents);

   // get nominal signal and background parameters

   Double_t signal_nominal = (1-0.0455)*gnorm*nevents/actual_integ;  // area outside 2sigma
   Double_t back_nominal = PMMks->Integral(mean-2*sig,mean+2*sig)*nevents/actual_integ - signal_nominal;
   Double_t StoB_nominal = signal_nominal/back_nominal;
   printf ("signal_nominal=%f, back_nominal=%f,StoB_nominal=%f\n",signal_nominal,back_nominal,StoB_nominal);

   Double_t ymin = 0;
   Double_t ymax = 300;
   // PMMks->GetYaxis()->SetRangeUser(ymin,ymax);
   // PMMks->DrawCopy();

   //
   TCanvas *c1 = new TCanvas("c1","c1 gen phithe",120,10,700,700);
   c1->SetBorderMode(0);
   c1->SetFillColor(0);
   c1->Divide(1,2);

   //c1->SetGridx();
   //c1->SetGridy();
   //c1->SetLogy();
   c1->cd(1);

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

   for (j=0;j<nevents;j++) {
      Double_t rMMks = PMMks->GetRandom();
      MMks->Fill(rMMks);
      }

   // get parameters of histogram

   Int_t ndim = MMks->GetNbinsX();
   Double_t xlo = MMks->GetBinLowEdge(1);
   Double_t width = MMks->GetBinWidth(1);
   printf ("gen_phithe: ndim=%d, xlo=%f, width=%f\n",ndim,xlo,width);

   // set starting parameters 
   Double_t mu=5;            // unis of 1/x
   Double_t sigma=0.04;      // units of x  
   Double_t gain=20;         // arbitrary
   Double_t xshift=1.47;     // x translation
   Double_t mean=1.538;       // gauss mean  
   Double_t sig=0.005;       // gauss sig  
   Double_t gnorm=40*sqrt(2*3.14159)*sig;       // gauss normalization 

   // get normalization in window
   Double_t window_norm=1;
   PMMks->SetParameters(mu,sigma,1,xshift,1,mean,sig);
   window_norm = PMMks->Integral(xmin,xmax)/2;
   printf (" window_norm=%f\n",window_norm);
   PMMks->SetParameters(mu,sigma,gain,xshift,gnorm,mean,sig);

   MMks->Fit("PMMks","","",xmin,xmax);
   Double_t mu_theta = PMMks->GetParameter(0);
   Double_t sigma_theta = PMMks->GetParameter(1);
   Double_t gain_theta = PMMks->GetParameter(2);
   Double_t xshift_theta = PMMks->GetParameter(3);
   Double_t gnorm_theta = fabs(PMMks->GetParameter(4));
   Double_t mean_theta = PMMks->GetParameter(5);
   Double_t sig_theta = fabs(PMMks->GetParameter(6));
   Double_t gnormerr_theta = PMMks->GetParError(4);
   Double_t gainerr_theta = PMMks->GetParError(2);
   Double_t splusb = PMMks->Integral(mean_theta-2*sig_theta,mean_theta+2*sig_theta)*nevents/PMMks->Integral(xmin,xmax);
   Double_t int_the = PMMks->Integral(xmin,xmax);
   printf ("gnorm_theta=%f, gain_theta*window_norm=%f, g+g/binsize=%f, int_the=%f, int_the/binsize=%f\n",gnorm_theta,gain_theta*window_norm,(gnorm_theta+gain_theta*window_norm)/binsize,int_the,int_the/binsize);

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

   // compute log likelihoods including signal
   Double_t logLBS, logLB;

   for(j=0;j<ndim;j++) {
      Double_t xbin = MMks->GetBinCenter(j+1);
      Double_t content = MMks->GetBinContent(j+1);
      Double_t func_theta = PMMks->Eval(xbin);
      logLBS = logLBS -func_theta + content*log(func_theta);
      printf ("j=%d, xbin=%f, content=%f, func=%f\n",j,xbin,content,func_theta);
   }
   printf ("logLBS=%f\n",logLBS);

   //
   // TCanvas *c2 = new TCanvas("c2","c2 gen phithe",200,10,700,700);
   c1->cd(2);
   c1_2->SetBorderMode(0);
   c1_2->SetFillColor(0);
   //c2->SetGridx();
   //c2->SetGridy();
   //c2->SetLogy();

   PMMks->FixParameter(4,0);
   PMMks->FixParameter(5,mean);
   PMMks->FixParameter(6,sig);
   MMks->Fit("PMMks","","",xmin,xmax); 
   Double_t mu_back = PMMks->GetParameter(0);
   Double_t sigma_back = PMMks->GetParameter(1);
   Double_t gain_back = PMMks->GetParameter(2);
   Double_t xshift_back = PMMks->GetParameter(3);
   Double_t gnorm_back = PMMks->GetParameter(4);
   Double_t mean_back = PMMks->GetParameter(5);
   Double_t sig_back = PMMks->GetParameter(6);
   Double_t b = PMMks->Integral(mean_theta-2*sig_theta,mean_theta+2*sig_theta)*nevents/PMMks->Integral(xmin,xmax);
   printf ("Nevents=%d, b =%f\n",nevents,b);

   // compute log likelihood for background only

   for(j=0;j<ndim;j++) {
      Double_t xbin = MMks->GetBinCenter(j+1);
      Double_t content = MMks->GetBinContent(j+1);
      Double_t func_back = PMMks->Eval(xbin);
      logLB = logLB -func_back + content*log(func_back);
      printf ("j=%d, xbin=%f, content=%f, func=%f\n",j,xbin,content,func_back);
   }
   printf ("logLB=%f\n",logLB);

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

   if (logLBS-logLB > 0) {
      Double_t Signif_log = sqrt(2*(logLBS-logLB));
      }
   else {
      Double_t Signif_log = 0;
      }

   printf ("Signif_log=%f, logLB =%f, logLBS=%f\n",Signif_log,logLB,logLBS);
   sprintf(string,"ln(LB)=%.1f\n",logLB);

   Double_t signal = (1-0.0455)*gnorm_theta/binsize;  // area outside 2sigma
   Double_t sigerr = (1-0.0455)*gnormerr_theta/binsize;
   Double_t back = PMMks->Integral(mean_theta-2*sig_theta,mean_theta+2*sig_theta)/binsize;
   Double_t backerr = back*gainerr_theta/gain_theta;
   if (backerr < sqrt(back)) backerr = sqrt(back);
   Double_t Signif = signal/backerr;
   printf ("signal=%f, sigerr=%f, sig_theta=%f, back=%f, backerr=%f,Signif=%f\n",signal,sigerr,sig_theta,back,backerr,Signif);

   c1->cd(1);
   t1 = new TLatex(0.45,0.85,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"ln(LB+S)=%.1f\n",logLBS);
   t1 = new TLatex(0.45,0.80,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"Log Sig=%.1f\n",Signif_log);
   t1 = new TLatex(0.45,0.75,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"Nevents=%d\n",nevents);
   t1 = new TLatex(0.45,0.70,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"Gen signal=%.3f\n",signal_nominal);
   t1 = new TLatex(0.45,0.65,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"Gen back=%.3f\n",back_nominal);
   t1 = new TLatex(0.45,0.60,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"Gen S/B=%.3f\n",StoB_nominal);
   t1 = new TLatex(0.45,0.55,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();

   // estimate the statistical error on the background under the peak

   Double_t s = (splusb-b)*(1-0.0455);
   Double_t signif_naive = s/sqrt(b);
   printf("splusb=%f, signal=%f, back=%f, Signif=%f, s=%f, b=%f, signif_naive=%f\n",splusb,signal,back,Signif,s,b,signif_naive);
   printf ("binsize=%f\n",binsize);

   sprintf(string,"Gen S/sqrt(B)=%.3f\n",signal_nominal/sqrt(back_nominal));
   t1 = new TLatex(0.45,0.50,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"Phi_norm=%.3f\n",phi_norm);
   t1 = new TLatex(0.65,0.55,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"Theta_norm=%.3f\n",theta_norm);
   t1 = new TLatex(0.65,0.50,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"Fit Signal=%.1f#pm%.1f\n",signal,sigerr);
   t1 = new TLatex(0.65,0.45,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"Fit Back=%.1f#pm%.1f\n",back,backerr);
   t1 = new TLatex(0.65,0.40,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"Fit S/B=%.3f\n",signal/back);
   t1 = new TLatex(0.65,0.35,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"S/sqrt(B)=%.3f\n",Signif);
   t1 = new TLatex(0.65,0.30,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();

   sprintf(filename,"gen_phithe_%d_%d.eps",100*phi_norm,100*theta_norm);
   printf ("file selection =%s\n",filename);
   c1->SaveAs(filename);
   // sprintf(filename,"gen_phithe_c2.eps",dataname);
   // c2->SaveAs(filename);
      

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

   Double_t mu=par[0];
   Double_t sigma=par[1];
   Double_t gain=par[2];
   Double_t xshift=par[3];
   Double_t norm=par[4];
   Double_t mean=par[5];
   Double_t sig=par[6];

   Double_t x1=x[0]-xshift;
   Double_t x2=x[0];         // use absolute x for gaussian

   char string[256];

   Double_t arg = (mu*sigma*sigma-x1)/(sqrt(2)*sigma);
   Double_t amplitude = (mu/2)* exp(-mu*x1 +(mu*sigma)*(mu*sigma)/2) *  TMath::Erfc(arg);
   Double_t gaus = norm*exp(-(x2-mean)*(x2-mean)/(2*sig*sig))/(sqrt(2*3.14159)*sig);

   amplitude = gain*amplitude + gaus;
   
   /*sprintf (string,"amplitude=%f mu=%f sigma=%f arg=%f\n",amplitude,mu,sigma,arg);
   printf ("string=%s",string);*/

   return amplitude;
}