void TurnOnFit(TH1 *h1, TH1 *h2, TH1 *h) {
  // Calculate the efficiency bin by bin as the ratio of two histograms.
  // Errors on the efficiency are calculated properly.

  int nbins = h->GetNbinsX();

  // Check histogram compatibility
  if (nbins != h1->GetNbinsX() || nbins != h2->GetNbinsX()) {
    std::cout << "Attempt to divide histograms with different number of bins" << std::endl;
    return;
  }

  // Issue a Warning if histogram limits are different
  if (h->GetXaxis()->GetXmin() != h1->GetXaxis()->GetXmin() ||
      h->GetXaxis()->GetXmax() != h1->GetXaxis()->GetXmax() ){
    std::cout << "Attempt to divide histograms with different axis limits" << std::endl;
  }

  if (h->GetXaxis()->GetXmin() != h2->GetXaxis()->GetXmin() ||
      h->GetXaxis()->GetXmax() != h2->GetXaxis()->GetXmax() ){
    std::cout << "Attempt to divide histograms with different axis limits" << std::endl;
  }

  //if (h->GetDimension() < 2) nbinsy = -1;

  // Loop on bins (including underflows/overflows)
  int bin;
  double b1,b2,effi,erro;
  for (bin=0;bin<=nbins+1;bin++) {
    
    b1  = h1->GetBinContent(bin);
    b2  = h2->GetBinContent(bin);
    //std::cout << "b1: " << b1 << endl;
    //std::cout << "b2: " << b2 << endl;

    if ( b2 ) {
      effi = b1/b2;
      erro = TMath::Sqrt((1-b1/b2)/b2*(b1/b2));
      //std::cout << "effi: " << effi << endl;
      //std::cout << "erro: " << erro << endl;
    } else {
      effi = 0;
      erro = 0;
    }
    h->SetBinContent(bin,effi);
    h->SetBinError(bin,erro);
  }
}

Double_t turnon_func(Double_t *x, Double_t *par)
{
  double halfpoint = par[0];
  double slope = par[1];
  double plateau = par[2];
  
  //double offset = par[3];
  //double plateau = 1.0;
  double offset = 0;

  double pt = TMath::Max(x[0],0.000001);

   double arg = 0;
   //cout << pt <<", "<< halfpoint <<", " << slope <<endl;
   arg = (pt - halfpoint)/(TMath::Sqrt(pt)*slope);
   double fitval = offset+0.5*plateau*(1+TMath::Erf(arg));
   return fitval;
}

Double_t turnon_func2(Double_t *x, Double_t *par)
{
  return turnon_func(x,par) + turnon_func(x,&par[3]);
}