Ok great, that looks awesome, but how do we know how strong our signal excess is?
For this, we can do a signal + background fit. The signal looks fairly Gaussian, so let's try that
TF1 *fit_sig_bkg = new TF1("fit_sig_bkg", "[0] + exp([2]*x+[1]) + [3]*TMath::Gaus(x,[4],[5])", 100, 160);
In order to help the fit converge, now that it has quite a lot of parameters, we can initialise the parameters to values that are about right
fit_sig_bkg->SetParameter(0, fit->GetParameter(0));
fit_sig_bkg->SetParameter(1, fit->GetParameter(1));
fit_sig_bkg->SetParameter(2, fit->GetParameter(2));
fit_sig_bkg->SetParameter(4, 125);
fit_sig_bkg->SetParameter(5, 3);
Now we can do the fit and draw it on our plot
pad1 -> cd();
fit_sig_bkg -> SetLineColor(kBlack);
h_data -> Fit(fit_sig_bkg);
fLegend -> AddEntry(fit_sig_bkg, "Fit signal+bkg", "l");
fLegend -> Draw(); // redraw the legend
We can pull the fit parameters using the GetParameter() function
cout << "Fitted signal normalisation: " << fit_sig_bkg->GetParameter(3) << " +- " << fit_sig_bkg->GetParError(3) << endl;
cout << "Fitted signal mean: " << fit_sig_bkg->GetParameter(4) << " +- " << fit_sig_bkg->GetParError(4) << endl;
cout << "Fitted signal width: " << fit_sig_bkg->GetParameter(5) << " +- " << fit_sig_bkg->GetParError(5) << endl;
The integral over all x of a gaussian is (exercise for student :-) ) sqrt(pi). For a Gaussian of width s
and normalisation n
it is therefore just s * n * sqrt(pi)
. So, we can calculate the number of signal events:
float fittedSignalEvents = fit_sig_bkg->GetParameter(3) * fit_sig_bkg->GetParameter(5) * TMath::Sqrt(TMath::Pi());
Now we can compare this to the number of signal events we expect from simulation. The signal strength mu is defined as the number of events we extract above our background prediction divided by the number of signal events predicted by our model. It should be close to 1 if we are seeing the same thing as our prediction!
float mu = fittedSignalEvents / hSig -> Integral();
cout << "Signal strength: " << mu << endl;
To round things off, let's plot the fitted signal in the ratio plot so we can see it by eye
pad2 -> cd();
TF1 * fit_sig = new TF1("fit_sig","gaus(0)",100, 160);
fit_sig->SetParameter(0, fit_sig_bkg->GetParameter(3));
fit_sig->SetParameter(1, fit_sig_bkg->GetParameter(4));
fit_sig->SetParameter(2, fit_sig_bkg->GetParameter(5));
fit_sig->SetLineColor(kBlack);
fit_sig->Draw("same");
There we go - we've discovered the Higgs! Of course, the actual analysis was a bit more complicated than this - in the next exercise we explore part of how the invariant mass distributions are created.