const bool scaleNumberOfBins = false;
const int N = 10; // number of widths
const double Wmax = 30.0; // maximum width
const RooFFTConvPdf::BufStrat strat = RooFFTConvPdf::Extend; // or Flat or Mirror

TCanvas* GaussianBias() {

  // Loop over range sizes
  double *rans = new double[N];
  double *vals = new double[N];
  double *errs = new double[N];

  double mult = Wmax / (2.0 * N);
  for (int i = 1; i <= N; i++) {

    // Generation and Fit models
    double range = i*mult;
    RooRealVar mass("mass","mass",0.0,-1*range,range);
    if (scaleNumberOfBins)
      mass.setBins(200 * i);
    else
      mass.setBins(2000);
    RooRealVar mean("mean","mean",0.0,-0.2,0.2);
    RooRealVar width("width","width",1.0,1e-3,2.0);
    width.setConstant(true);
    RooGaussian truth("truth","truth",mass,mean,width);
    RooConstVar bias("bias","bias",0.0);
    RooRealVar res("res","res",1.0,1e-3,2.0);
    res.setConstant(true);
    RooGaussian smear("smear","smear",mass,bias,res);
    RooFormulaVar totalWidth("totalWidth","sqrt(width*width + res*res)",RooArgList(width,res));
    RooGaussian obs_ana("obs_ana","obs_ana",mass,mean,totalWidth);
    RooFFTConvPdf obs_num("obs_num","obs_num",mass,truth,smear);
    obs_num.setBufferStrategy(strat);

    // Generate data from analytical pdf
    RooDataSet* data = obs_ana.generate(mass,5000000);
    obs_num.fitTo(*data);

    // Save val
    rans[i-1] = 2*range;
    vals[i-1] = mean.getVal();
    errs[i-1] = mean.getError();
  }
  
  // Make graph
  TGraphErrors* graph = new TGraphErrors(N,rans,vals,0,errs);
  graph->SetName("bgraph");
  graph->SetTitle("bgraph");
  graph->SetLineColor(2);
  graph->SetLineWidth(4);
  graph->SetMarkerColor(4);
  graph->SetMarkerStyle(21);
  graph->GetXaxis()->SetTitle("Range size");
  graph->GetYaxis()->SetTitle("Bias in numerical fit");

  // Make and return canvas
  TCanvas* c = new TCanvas("canvas","canvas",800,800);
  c->cd();
  graph->Draw("AP");
  TF1 f("f","[0]*x",0,2*Wmax);
  graph->Fit(&f);
  c->Update();
  return c;
}