void confidenceBelt()
{
   // the model is a Voigt function. g=2, sigma=1+x/10
   // we scan x from 20 to 100 by steps of 1
   // for each case, we compute the low and high edges
   // We create a single graph for with the confidence belt
   Double_t x[290];
   Double_t y[290];
   TF1* f = new TF1("myvoigt","TMath::Voigt(x-[0],[1],[2])",-50,200);
   f->SetNpx(10000);
   Double_t q[2];
   Double_t probsum[2] = {0.15865,0.84135};
   int i=0;
   for(float x0 = 5; x0<150; x0+=1) {
     f->SetParameters(x0,1+x0/10.,2);
     f->SetRange(x0-5*(x0+3),x0+5*(x0+3)); // we need this to regularize the quantile algo
     f->GetQuantiles(2,q,probsum); // we use here (approximation of) central intervals
     y[i] = x0;
     x[i] = q[0];
     y[289-i] = x0;
     x[289-i] = q[1];
     i++;
   }
   TGraph* g = new TGraph(290,x,y);
   g->SetFillColor(kBlue);
   g->SetTitle("Confidence Belt");
   g->Draw("AF");
   g->GetXaxis()->SetRangeUser(20,100);

   // now, imagine we do a measurement at 45, what is the 68% confidence interval
   g->GetXaxis()->SetRangeUser(20,130);
   Float_t measurement = 45.;
   TLine* line = new TLine(measurement,0,measurement,160);
   line->Draw("same");
   TGraph* subg1 = new TGraph(145,x,y);
   TGraph* subg2 = new TGraph(145,x+145,y+145);
   std::cout << "For a measurement at " << measurement << " the confidence interval is [" << subg2->Eval(measurement) << ", " << subg1->Eval(measurement) << "]. " << std::endl; 
   
}