double Neff ( double *x, double *par ) {
  double phi = x[0];

  double c = par[0];
  double gc = par[1];
  double Nc0 = par[2];
  double Neff0 = par[3];

  double Nc = 0;

  Nc = Nc0 * (1. - TMath::Exp( - c * phi ) ) + gc * phi;

  double Neff = Neff0 - Nc;

  return fabs(Neff);
}


void deltaNeff( double Neff0 = 1.1e12  ) {
  double deltaNeff = 0;

  double Nc0 = 9.0e11;
  double c  = 1e-16;
  double gc = 1e-2;

  double phimin = 1e11;
  double phimax = 2e16;

  double phi = 0;

  double *x = φ

  double par[4] = { c, gc, Nc0, Neff0 };

  TF1 *f = new TF1("f",Neff,phimin,phimax,4);
  f->SetParameters( par );
  f->SetParNames( "c", "g_{c}", "N_{c,0}", "N_{eff,0}" );

  TCanvas *c1 = new TCanvas("c1","",800,600);
  c1->cd();

  f->Draw("l");

  f->GetYaxis()->SetTitle("|N_{eff}| [cm^{-3}]");
  f->GetXaxis()->SetTitle("#phi [n_{eq}/cm^{-2}]");

  gPad->SetLogx(1);
  gPad->SetLogy(1);

  f->SetTitle("|N_{eff}| & V_{depl} vs #phi");

  double e0 = 8.85e-12/100;
  double er = 11.9;
  double  q = 1.6e-19;
  double  w = 300e-4;

  double neffmin = f->GetMinimum(phimin,phimax);
  double neffmax = f->GetMaximum(phimin,phimax);

  neffmin = 1.0e11;
  std::cout << "Neff min = " << neffmin << "\n";
  std::cout << "Neff max = " << neffmax << "\n";

  double vmin = neffmin*q*w*w/2/e0/er ;
  double vmax = neffmax*q*w*w/2/e0/er;
 


   //draw an axis on the right side
   TGaxis *axis = new TGaxis(phimax,neffmin,phimax,neffmax,
         vmin,vmax,510,"+G");
   axis->SetLineColor(kRed);
   axis->SetLabelColor(kRed);
   axis->SetTextColor(kRed);
   axis->SetTitle("V_{depl} [V]");
   axis->SetLabelOffset(0.04);
   axis->Draw();

   gPad->SetGrid(2,2);

   gPad->Update();
   gPad->Modified();

   const int nmeas = 11;
   double xbins[nmeas] = {1e12,1e13,2e13,5e14,1e15,2e15,3e15,4e15,5e15,,7e15,1e16};
   double ybins[nmeas] = {
     f->Eval(1e12)*(1.+gRandom->Gaus()*0.05),
     f->Eval(1e13)*(1.+gRandom->Gaus()*0.05),
     f->Eval(2e13)*(1.+gRandom->Gaus()*0.05),
     f->Eval(5e14)*(1.+gRandom->Gaus()*0.05),
     f->Eval(1e15)*(1.+gRandom->Gaus()*0.05),
     f->Eval(2e15)*(1.+gRandom->Gaus()*0.05),
     f->Eval(3e15)*(1.+gRandom->Gaus()*0.05),
     f->Eval(4e15)*(1.+gRandom->Gaus()*0.05),
     f->Eval(5e15)*(1.+gRandom->Gaus()*0.05),
     f->Eval(7e15)*(1.+gRandom->Gaus()*0.05),
     f->Eval(1e16)*(1.+gRandom->Gaus()*0.05)
   };
   double ex[nmeas] = {xbins[0]*.15,xbins[1]*.15,xbins[2]*.15};
   double ey[nmeas] = {ybins[0]*.15,ybins[1]*.15,ybins[2]*.15};
   TGraphErrors *Irr = new TGraphErrors(nmeas,xbins,ybins,ex,ey);

   Irr->SetName("Irr");
   Irr->SetTitle("Irradiation data");

   gPad->Update();
   gPad->Modified();

   Irr->SetMarkerStyle(20);
   Irr->Draw("sameP");

   Irr->Fit(f);
//
   gPad->Update();
   gPad->Modified();

   Irr->GetYaxis()->SetTitle("|N_{eff}| [cm^{-3}]");
   Irr->GetXaxis()->SetTitle("#phi [n_{eq}/cm^{-2}]");
//
   gPad->Update();
   gPad->Modified();

   TPaveStats* ps = (TPaveStats*)Irr->GetListOfFunctions()->FindObject("stats");
   ps->SetX1NDC(0.15);
   ps->SetY1NDC(0.57);
   ps->SetX2NDC(0.51);
   ps->SetY2NDC(0.85);

   gPad->Update();
   gPad->Modified();

   c1->SaveAs("deltaNeff.png");

  std::cout << "NDF = " << f->GetNDF() << "\n\n\n\n\n";

}