Problem fitting Toy mc

Newbie
1

Hello,
I am trying to do a toy mc with fixed bins, where I have a 5 parameter function, which I use to generate data and to fit the data that it generates. For some reason when I run this script, the p-value is 0. I don’t understand why this fit is so bad. Can someone explain? Is there a way to make it work?

TCanvas* c1 = new TCanvas("c","BPRE",10,10,600,600);
c1->Divide(1,1,0.005,0.005);
c1->SetTickx();
c1->SetTicky();
c1->SetTitle("");
c1->SetLineWidth(3);
c1->SetBottomMargin(0.1);
c1->SetTopMargin(0.05);
c1->SetRightMargin(0.01);
c1->SetFillColor(0);
TPad* pad1 = new TPad("pad1","pad1",0,0.3,1,0.97);pad1->SetBottomMargin(0);
pad1->SetLeftMargin(0.13);
pad1->SetRightMargin(0.04);
pad1->SetTopMargin(0.02);
pad1->Draw();
pad1->cd();
pad1->SetLogy(1);
pad1->SetLogx(1);
TH1F *h=pad1->DrawFrame(0.387,.9,10.1,100000000.0-100000.0);
TAxis *ay = h->GetYaxis();
ay->SetTitle("Events/TeV");
TF1 parabola_5("para_five","([0]*((1-(x/13.0))**[1])*((x/13.0)**([2]+([3]*TMath::Log((x/13.0)))+([4]*(TMath::Log((x/13.0))^2)))))",0.387,10.1);
parabola_5.SetParameter(0,4.03063e+06);
parabola_5.SetParameter(1,1.67322e+01);
parabola_5.SetParameter(2,7.27028e+00);
parabola_5.SetParameter(3,3.53873e+00);
parabola_5.SetParameter(4,4.62872e-01);

double edges[132]={0.387,0.405,0.424,0.443,0.462,0.482,0.502,0.523,0.544,0.566,0.588,0.611,0.634,0.657,0.681,0.705,0.73,0.755,0.781,0.807,0.834,0.861,0.889,0.917,0.946,0.976,1.006,1.037,1.068,1.1,1.133,1.166,1.2,1.234,1.269,1.305,1.341,1.378,1.416,1.454,1.493,1.533,1.573,1.614,1.656,1.698,1.741,1.785,1.83,1.875,1.921,1.968,2.016,2.065,2.114,2.164,2.215,2.267,2.32,2.374,2.429,2.485,2.542,2.6,2.659,2.719,2.78,2.842,2.905,2.969,3.034,3.1,3.167,3.235,3.305,3.376,3.448,3.521,3.596,3.672,3.749,3.827,3.907,3.988,4.07,4.154,4.239,4.326,4.414,4.504,4.595,4.688,4.782,4.878,4.975,5.074,5.175,5.277,5.381,5.487,5.595,5.705,5.817,5.931,6.047,6.165,6.285,6.407,6.531,6.658,6.787,6.918,7.052,7.188,7.326,7.467,7.61,7.756,7.904,8.055,8.208,8.364,8.523,8.685,8.85,9.019,9.191,9.366,9.544,9.726,9.911,10.1};
auto pfivehisto = new TH1F("5p background","Signal and Background;xvals;Events / TeV",131,edges);

pfivehisto->SetMarkerStyle(8);
pfivehisto->SetLineColor(kRed);
pfivehisto->SetTitle("Signal and Background");
pfivehisto->Scale(1.0/92);
gStyle->SetOptFit(1111);
TRandom3 *r = new TRandom3();
r->SetSeed(time(0));
for (int i=1; i<=4.65903e+06;++i){pfivehisto->Fill(parabola_5.GetRandom(0.387,10.1));}
int npar = parabola_5.GetNpar();
cout << "npar=" << npar << endl;
for (int i=1; i<100; i++){pfivehisto->Fit(&parabola_5,"SL","SAMES");}

_ROOT Version: 6.23/01

2

Hi,

The problem is due to the variable width histogram.
You have two possibility:

  1. Scale the histogram by dividing the content by the bin width before fitting. The drawback of doing this is that errors in the maximum likelihood fit (option L) will be approximated (you need to use LW)
TRandom3 *r = new TRandom3();
r->SetSeed(0); // why time (0) ? you need to use 0 for different seed every time
parabola_5->SetNpx(1000); // important to get good sampling from TF1
for (int i=1; i<=4.65903e+06;++i){pfivehisto->Fill(parabola_5.GetRandom(0.387,10.1));}
pfivehisto->Scale(1.,"width");
pfivehisto->Fit(&parabola_5,"SLW","SAMES");
  1. better solution, use fit option WIDTH. It requires to re-drawn the histogram afterwards. (This maybe should be fixed on the ROOT side and automatically scale the histogram)
TRandom3 *r = new TRandom3();
r->SetSeed(0); // why time (0) ? you need to use 0 for different seed every time
parabola_5->SetNpx(1000); // important to get good sampling from TF1
for (int i=1; i<=4.65903e+06;++i){pfivehisto->Fill(parabola_5.GetRandom(0.387,10.1));}
pfivehisto->Fit(&parabola_5,"SL WIDTH","");
pfivehisto->Scale(1.,"width");
pfivehisto->Draw(); 
parabola_5.Draw("SAME");

Lorenzo

1 Like
3

Ah ok, thank you very much!

4

Now I have a new problem, I plot the residuals, but now for some reason, the negative residuals do not show on the graph. I assume it has something to do with the variable bin width because when I used the fixed bin width I did not encounter this problem.

TCanvas* c1 = new TCanvas("c","BPRE",10,10,600,600);
c1->Divide(1,1,0.005,0.005);
c1->SetTickx();
c1->SetTicky();
c1->SetTitle("");
c1->SetLineWidth(3);
c1->SetBottomMargin(0.1);
c1->SetTopMargin(0.05);
c1->SetRightMargin(0.01);
c1->SetFillColor(0);
TPad* pad1 = new TPad("pad1","pad1",0,0.3,1,0.97);pad1->SetBottomMargin(0);
pad1->SetLeftMargin(0.13);
pad1->SetRightMargin(0.04);
pad1->SetTopMargin(0.02);
pad1->Draw();
pad1->cd();
pad1->SetLogy(1);
pad1->SetLogx(1);
TH1F *h=pad1->DrawFrame(0.387,.00009,10.1,100000000.0-100000.0);
TAxis *ay = h->GetYaxis();
ay->SetTitle("Events/TeV");
TF1 parabola_5("para_five","([0]*((1-(x/13.0))**[1])*((x/13.0)**([2]+([3]*TMath::Log((x/13.0)))+([4]*(TMath::Log((x/13.0))^2)))))",0.387,10.1);
parabola_5.FixParameter(0,1.51746e+08);
parabola_5.FixParameter(1,16.6014);
parabola_5.FixParameter(2,7.1695);
parabola_5.FixParameter(3,3.51143);
parabola_5.FixParameter(4,0.460357);

double edges[132]={0.387,0.405,0.424,0.443,0.462,0.482,0.502,0.523,0.544,0.566,0.588,0.611,0.634,0.657,0.681,0.705,0.73,0.755,0.781,0.807,0.834,0.861,0.889,0.917,0.946,0.976,1.006,1.037,1.068,1.1,1.133,1.166,1.2,1.234,1.269,1.305,1.341,1.378,1.416,1.454,1.493,1.533,1.573,1.614,1.656,1.698,1.741,1.785,1.83,1.875,1.921,1.968,2.016,2.065,2.114,2.164,2.215,2.267,2.32,2.374,2.429,2.485,2.542,2.6,2.659,2.719,2.78,2.842,2.905,2.969,3.034,3.1,3.167,3.235,3.305,3.376,3.448,3.521,3.596,3.672,3.749,3.827,3.907,3.988,4.07,4.154,4.239,4.326,4.414,4.504,4.595,4.688,4.782,4.878,4.975,5.074,5.175,5.277,5.381,5.487,5.595,5.705,5.817,5.931,6.047,6.165,6.285,6.407,6.531,6.658,6.787,6.918,7.052,7.188,7.326,7.467,7.61,7.756,7.904,8.055,8.208,8.364,8.523,8.685,8.85,9.019,9.191,9.366,9.544,9.726,9.911,10.1};
auto pfivehisto = new TH1F("5p background","Signal and Background;xvals;Events / TeV",131,edges);

pfivehisto->SetMarkerStyle(8);
pfivehisto->SetLineColor(kRed);
pfivehisto->SetTitle("Signal and Background");
pfivehisto->Scale(1.0/92);
gStyle->SetOptFit(1111);
TRandom3 *r = new TRandom3();
r->SetSeed(0); //time(0)
parabola_5.SetNpx(1000); // important to get good sampling from TF1
for (int i=1; i<=4.65903e+06;++i){pfivehisto->Fill(parabola_5.GetRandom(0.387,10.1));}
int npar = parabola_5.GetNpar();
cout << "npar=" << npar << endl;
pfivehisto->Fit(&parabola_5,"SL WIDTH","");
pfivehisto->Scale(1.,"width");
pfivehisto->SetAxisRange(.000001,100000000,"Y");
pfivehisto->Draw("SAME");
parabola_5.Draw("SAME");
TF1* g = new TF1("gaus","[0]*TMath::Gaus(x,[1],[2])",0.721545,1.03845);
double a=0.1, b=0.88, c=0.031691;
g->SetParameter(0,a);
g->SetParameter(1,b);
g->SetParameter(2,c);
g->SetLineColor(kBlue);
//ay->SetRangeUser(.00001,100000000);
g->Draw("Same");
for (int i=1; i<=(0.5*1964);++i) pfivehisto->Fill(g->GetRandom(0.721545,1.03845));
for (int i=1; i<10; i++){auto fitr = pfivehisto->Fit(&parabola_5,"SL","SAMES");fitr->Print();}
pfivehisto->Draw("PE SAME");
auto leg = new TLegend(.1,.7,.5,.9,"");
leg->SetFillColor(0);
leg->SetTextSize(0.02);
leg->AddEntry(pfivehisto,"Background Points and Gaussian Points");
leg->Draw();
leg->AddEntry(&parabola_5,"Background Fit (Fixed)");
leg->Draw();
c1->Update();
c1->cd();
TPad* pad2 = new TPad("pad2","pad2",0,0.01,1,0.3);
pad2->SetLeftMargin(0.13);
pad2->SetRightMargin(0.04);
pad2->SetTopMargin(0);
pad2->SetBottomMargin(0.31);
pad2->Draw();
pad2->cd();
TH1F *h1=pad2->DrawFrame(.4,-5,10.0,5);
TAxis *ax = h1->GetXaxis();
ay = h1->GetYaxis();
ay=h1->GetYaxis();
ay->SetRangeUser(-5,5);
ay->SetLabelFont(42);
ay->SetLabelSize(0.10);
ay->SetTitleSize(0.3);
ay->SetNdivisions(505);
ax=h1->GetXaxis();
ax->SetTitle("m_{jjl} (TeV)");
ax->SetTitleOffset(1.18);
ay->SetTitleOffset(0.8);
ax->SetLabelFont(42);
ay->SetLabelFont(42);
ax->SetLabelSize(0.12);
ax->SetTitleSize(0.14);
ay->SetLabelSize(0.12);
ay->SetTitleSize(0.2);
ax->Draw("same");
ay->Draw("same");
TGraphErrors DiffBG=TGraphErrors();
DiffBG.SetMarkerColor(2);
DiffBG.SetMarkerStyle( 20 );
DiffBG.SetMarkerSize( 0.7 );
TH1 *Sig= (TH1*) pfivehisto->Clone();
Sig->SetFillColor(2);
TH1D res=TH1D("Res","Res",100,-10,10);
int n=0;
for (int i=1; i<pfivehisto->GetNbinsX(); i++){Sig->SetBinContent(i,0);Sig->SetBinError(i,0);if(pfivehisto->GetBinCenter(i)>0.387&&pfivehisto->GetBinCenter(i)<10.0 && pfivehisto->GetBinContent(i)>0){double center=pfivehisto->GetBinCenter(i);double x=pfivehisto->GetBinCenter(i);double D =pfivehisto->GetBinContent(i);double Derr = pfivehisto->GetBinError(i);double B =parabola_5.Eval(center);double Berr = sqrt(B);double frac=0.0;double fracErr = 0.0;if (B!=0){frac = (D-B)/Derr;cout << frac<<endl;}res.Fill(frac);Sig->SetBinContent(i,frac);Sig->SetBinError(i,fracErr);DiffBG.SetPoint(n,x,frac);n++;}}
TLine ar5=TLine(c1->XtoPad(.4),0,c1->XtoPad(10),0);
ar5.SetLineWidth(2);
ar5.SetLineStyle(2);
ar5.Draw("same");
TLatex l5=TLatex();
l5.SetTextAngle(90);
l5.SetTextSize(0.10);
l5.SetNDC();
l5.SetTextColor(1);
l5.DrawLatex(0.03,0.4,"D_{i} - F_{i} / #Delta D_{i}");
Sig->Draw("same histo ][");
pad2->SetLogx(1);
5

Ok, I solved this problem, I just included the line

Sig->Scale(1.);

after the line

Sig->Draw("same histo ][");

and that solved the problem. Not sure why, but it works nonetheless.

6

Hmmm, I wonder how

can change anything. It just multiplies your TH1 by 1.

7

Well, it does, you can try running it with and without that line, and you will see it makes a difference. At least it is different when I run it.