/*******************************************************************************
  Subtracting background from data not taking into account overflows and underflows
*******************************************************************************/
void subtrBgdFromData(TH1D *a, TH1D *b, float &nbgdThr, float nbgdThr_err, float &nsigThr, float nsigThr_err, int isol_th, bool withOverflows)
{
  TH1D *data  = (TH1D *) a->Clone();
  TH1D *templ = (TH1D *) b->Clone();

  //  verifying axis limits
  /* 
    std::cout << "axis limits" << data->GetXaxis()->GetXmax() << " " << templ->GetXaxis()->GetXmax() << std::endl;
    std::cout << "axis limits" << data->GetXaxis()->GetXmin() << " " << templ->GetXaxis()->GetXmin() << std::endl;
  */ 

  int bin0 = 0;
  if( withOverflows ) bin0 = 0;
  // Find bin for threshold value
  int bin1=1;
  for(;bin1<20;bin1++) if(floor(data->GetBinLowEdge(bin1)*100+0.5)>=isol_th) break;
  --bin1;

  double maxX = data->GetXaxis()->GetXmax();
  int    bin2 = data->FindBin(maxX)-1;
  if( withOverflows ) bin2 =  data->FindBin(maxX);

  printf("bin1 = %4d, bin2=%4d \n",bin1, bin2);

  float nData_L = data->Integral(bin0, bin1);
  float nData_L_err = sqrt(nData_L);
  float nData_R = data->Integral(bin1+1, bin2);
  float nData_R_err = sqrt(nData_R);

  printf("nData_L = %5.2f+-%5.2f, nData_R = %5.2f+-%5.2f \n",nData_L, nData_L_err, nData_R, nData_R_err);

  //introduce template normalisation below and above threshold region
  float nTempl_L = templ->Integral(bin0, bin1);
  float nTempl_R = templ->Integral(bin1+1, bin2);
  // to avoid numerical instability
  if(  nTempl_L == 0 )  nTempl_L = 1.; 
  if(  nTempl_R == 0 )  nTempl_R = 1.;
  float nTempl_L_err = sqrt(nTempl_L);
  float nTempl_R_err = sqrt(nTempl_R);

  // introduce ratio of tempalte integral below and above threshold
  // allows to uncorrelate errors in error propagation 
  float beta = nTempl_L/nTempl_R;
  float beta_err = beta * sqrt( 1./nTempl_L + 1./nTempl_R);

  printf("nTempl_L = %5.2f+-%5.2f, nTempl_R = %5.2f+-%5.2f \n",nTempl_L, nTempl_L_err, nTempl_R, nTempl_R_err);
  printf("beta = %5.2f+-%5.2f \n",beta, beta_err);

  // assume right to the threshold only background in the data
  float nBgd_R = data->Integral(bin1+1,bin2);
  float nBgd_R_err = sqrt(nBgd_R);
  // calculate bgd left to the threshold using beta scaling factor
  float nBgd_L = beta * nBgd_R;
  float nBgd_L_err = sqrt(   (beta_err * nBgd_R)* (beta_err * nBgd_R) 
                           + (beta * nData_R_err)*(beta * nData_R_err));

  printf("nBgd_L = %5.2f+-%5.2f, nBgd_R = %5.2f+-%5.2f \n",nBgd_L, nBgd_L_err, nBgd_R, nBgd_R_err);

  //now calculate signal
  //assume right to the threshold only  background in the data
  float nSig_L = nData_L - beta * nBgd_R;
  float nSig_L_err = sqrt( nData_L_err* nData_L_err + nBgd_L_err*nBgd_L_err );

  float nSig = nSig_L;
  float nSig_err = nSig_L_err;

  float nBgd = (1 + beta) * nData_R;
  float nBgd_err = sqrt(   (1 + beta)*(1+beta) * nData_R_err* nData_R_err
			 + beta_err*beta_err   * nData_R*nData_R);


  printf("nSig_L = %5.2f+-%5.2f \n",nSig_L, nSig_L_err);


  nbgdThr = nBgd_L;
  nbgdThr = nBgd_L_err;

  nsigThr = nSig_L;
  nsigThr_err = nSig_L_err;

  delete data;
  delete templ;
}