void makedata(double *data, const int N){
  TRandom2 r,p,q;

  for(int ievent=0;ievent<N;ievent++){
    if (p.Uniform()<0.98){
      data[ievent] = r.Exp(3.);
    }else{
      data[ievent] = q.Gaus(2., 0.1);
    }
  }
}

void histdata(double *data, TH1F *histo, const int N){
  //fill histogream from an array of data pointts
  for (int j = 0; j<N; j++){
    histo->Fill(data[j]);
  }

}

double expgaus(double *x, double *par){
  //this is the pdf of signal+background (s+b) hypothesis scaled to
  //the data size N :
  //i.e. expgaus

  const double root2pi = 2.5066;
  double alpha = 0.02, binwidth = 0.1, N = 100000.;
  
  double m0 = 2.0, sigma = 0.1;
  double xx = x[0];
  double f = par[0]*exp(-xx/par[1]) + par[2]*N*alpha*binwidth*(1./(root2pi*sigma))*exp(-0.5*(xx-m0)*(xx-m0)/(sigma*sigma));
  return f;
} 


void fitexpbinned(){

  //example code for profile likelihood ratio calculation -sb

  const int N = 100000;
  double m_gg[N];

  //generate data
  makedata(m_gg, N);

  //bin the  data
  const int nbin = 100;
  TH1F *hdata = new TH1F("data", "exp distribution", nbin, 0., 10.);
  histdata(m_gg, hdata, N);

  //visualize binned data
  hdata->SetMarkerStyle(8);
  hdata->Draw();

  //fit binned data
  //A is the normalization and K is the mean of the exponential to fit 
  double A, K;
  double A_start = 100, A_step = 100., K_start = 1.0, K_step = 0.1;
  double ll_max = -10000000., K_max = -10000000., A_max = -10000000.;
  double n[nbin],x[nbin];

  //get the bin centers and bin contents in two arrays
  for (int ibin=1; ibin<nbin; ibin++){
    x[ibin] = hdata->GetBinCenter(ibin);
    n[ibin] = hdata->GetBinContent(ibin);
  } 

  //b only maximization 
  for (int iA=0; iA<100; iA++){
    A = A_start + iA*A_step;
    for (int iK=0; iK<100; iK++){
      K = K_start + iK*K_step;
      double ll = 0.;
      for (int ibin=1; ibin<nbin; ibin++){
        if (ibin<14||ibin>25){
	  double lambda = A*exp(-x[ibin]/K);
	  ll += -lambda + n[ibin]*log(lambda);
	}
        if (ll>ll_max) {
	  ll_max = ll;
          K_max  = K;
	  A_max  = A;
	}
      } 
      //cout<<"A= "<<A<<" K= "<<K<<" ll= "<<ll<<endl;
    }
  }

  //result of maximization
  cout<<"A_max= "<<A_max<<" K_max= "<<K_max<<" ll_max= "<<ll_max<<" (Side band only"<<endl;
     ll = 0.0;
       for (int ibin=1; ibin<nbin; ibin++){
	  double lambda = A_max*exp(-x[ibin]/K_max);
	  ll += -lambda + n[ibin]*log(lambda);
	}  

    cout<<"ll_max= "<<ll<<"  (Full range of data, background only)"<<endl;

  TF1 *fb = new TF1("bfit","[0]*exp(-x/[1])",0.,10.);
  fb->SetParameter(0,A_max);
  fb->SetParameter(1,K_max);
  fb->SetLineColor(4);   
  fb->Draw("same");
 

  //s + b maximization 
  //  for (int iA=0; iA<100; iA++){
  //    A = A_start + iA*A_step;
  //    for (int iK=0; iK<100; iK++){
  //     K = K_start + iK*K_step;
  double mu, mu_start = 0.0, mu_step = 0.02,mu_max = -1.0;
  double par[3];
  par[0] = A_max;
  par[1] = K_max;

  double mu = mu_start;

  TF1 *fsb = new TF1("sbfit",expgaus,0.,10.,3);
  fsb->SetParameter(0,A_max);
  fsb->SetParameter(1,K_max);
  fsb->SetParameter(2,mu);
  
  ll_max = -10000000.; //reset ll_max
  for (int im=0; im<100; im++){
    mu = mu_start + im*mu_step;
    par[2] = mu;
    double ll = 0.;
    for (int ibin=1; ibin<nbin; ibin++){
      //      double lambda = expgaus(x[ibin], par);
      fsb->SetParameter(2,mu);
      double lambda = fsb->Eval(x[ibin]);
      ll += -lambda + n[ibin]*log(lambda);
      if (ll>ll_max) {
	ll_max = ll;
        mu_max  = mu;
	A_max  = A;
      }
    } 
    //cout<<"A= "<<A<<" K= "<<K<<" ll= "<<ll<<endl;
  }
  cout<<"mu_max= "<<mu_max<<" "<<" ll_max= "<<ll_max<<endl;
  
  
  fsb->SetParameter(2,mu_max);
  fsb->SetLineColor(6);
  fsb->Draw("same");

}