Dear @jonas ,
Thank you very much for your response.
I’m using:
ROOT 6.30/06
RedHat Enterprise Linux 9
gcc version 11.4.1 20230605 (Red Hat 11.4.1-2) (GCC)
1-
Before I share my code I have to mention that I found the part of the code that causes the errors. In my code I try to calculate the Chi2/ndf and that’s where I see the errors.
Also, the Chi2/ndf calculated is very large. I know this because when I fit my data in root (not RooFit) to double-Gaussian I have a Chi2/ndf of about 5, but when I fit data with my Roofit code to a Double-Gaussian, I have a Chi2/ndf if around 10,000.
2-
My other problem is that when I try to get the values from the fit corresponding to different x values, the y values does not match the values that I see on plots or get from my histogram. I know this is because My model is a PDF and I have to use a kind of a scale factor. Probably if I could calculate the integral of my model I could use that to calculate the scale factor but I don’t know how. My model is a “RooAddPdf”
object. I don’t know how to calculate it’s integral.
Thank you very much for your help,
I removed unnecessary parts of the code to make it shorter.
Here is my code:
#include "RooFitter.h"
void RooFitter::FitResults()
{
float PDGMass = 0.0;
const float JPsiPDGMass = 3.097;
const float PsiPDGMass = 3.686;
const float BPDGMass = 5.279;
TString Branch;
TCut SelectionCut;
PDGMass = BPDGMass;
Branch = "BMass";
LeftCut = 5.16;
RightCut = 5.40;
NBins = (RightCut - LeftCut)/0.0025; // <---<< 2.5MeV/Bin (96 Bins)
SelectionCut = "((BMass>5.16 && BMass<5.20) || (BMass>5.36 && BMass<5.40))";
for (int i=1; i<=h_FitResult->GetNbinsX(); ++i) {
double BinContent = h_FitResult->GetBinContent(i);
double BinError = sqrt(BinContent); // Assuming Poisson uncertainties
h_FitResult->SetBinError(i, BinError);
}
RooRealVar Mass(Branch, Branch, LeftCut, RightCut);
RooArgSet Variables(Mass);
RooArgSet NormalizationSet(Mass);
RooDataHist DataHist("DataHist", "Data Histogram", Mass, Import(*h_FitResult));
RooDataSet* Data = new RooDataSet("Data", "Data", RooArgSet(Mass));
for (int i=0; i<h_FitResult->GetNbinsX(); ++i) {
Mass.setVal(h_FitResult->GetBinCenter(i+1));
for (int j=0; j<h_FitResult->GetBinContent(i+1); ++j) {
Data->add(RooArgSet(Mass));
}
}
RooDataSet* ReducedData = (RooDataSet*)Data->reduce(SelectionCut);
// Signal and Background Modeling
#if defined FixMass
RooRealVar Mean("Mean","Mean",PDGMass);
#else
RooRealVar Mean("Mean","Mean",PDGMass,PDGMass-.020,PDGMass+.020);
#endif
RooRealVar* Sigma1 = nullptr;
RooRealVar* Sigma2 = nullptr;
if (SignalType == "B") {
if (JorP == "J") {
Sigma1 = new RooRealVar("Sigma1","Sigma of the first Gaussian", 1.110e-2, 9.000e-3, 1.320e-2);
Sigma2 = new RooRealVar("Sigma2","Sigma of the second Gaussian", 2.500e-2, 2.000e-2, 3.000e-2);
}
else if (JorP == "P") {
Sigma1 = new RooRealVar("Sigma1","Sigma of the first Gaussian", 1.110e-2, 9.00e-3, 1.320e-2);
Sigma2 = new RooRealVar("Sigma2","Sigma of the second Gaussian", 2.500e-2, 2.000e-2, 3.000e-2);
}
}
else if (SignalType == "P") {
Sigma1 = new RooRealVar("Sigma1","Sigma of the first Gaussian", 1.110e-2, 9.00e-3, 1.320e-2);
Sigma2 = new RooRealVar("Sigma2","Sigma of the second Gaussian", 2.500e-2, 2.000e-2, 3.000e-2);
}
else if (SignalType == "J") {
Sigma1 = new RooRealVar("Sigma1","Sigma of the first Gaussian", 1.110e-2, 9.00e-3, 1.320e-2);
Sigma2 = new RooRealVar("Sigma2","Sigma of the second Gaussian", 2.500e-2, 2.000e-2, 3.000e-2);
}
// Double-Gaussian P.D.F.
RooGaussian Gaussian1("Gaussian1","First Gaussian P.D.F.", Mass, Mean, *Sigma1);
RooGaussian Gaussian2("Gaussian2","Second Gaussian P.D.F.", Mass, Mean, *Sigma2);
RooRealVar coef("coef", "coef", 0.5, 0, 1);
RooAddPdf Signal("Signal", "Double-Gaussian Signal", RooArgList(Gaussian1, Gaussian2), coef);
RooRealVar x("x", "x", LeftCut, RightCut);
int NPoints = 0;
if (SignalType == "B")
NPoints = 7;
else if (SignalType == "P" || SignalType == "J")
NPoints = 20;
TH1D* FNB = (TH1D*)h_FitResult->Clone("FNB");
TH1D* LNB = (TH1D*)h_FitResult->Clone("LNB");
FNB->GetXaxis()->SetRange(1, NPoints);
LNB->GetXaxis()->SetRange((NBins-NPoints+1), NBins);
TH1D* CombinedHist = (TH1D*)FNB->Clone("CombinedHist");
CombinedHist->Add(LNB);
RooDataHist CombinedData("CombinedData", "Combined Data from first and last N bins", x, RooFit::Import(*CombinedHist));
RooRealVar m("m", "m", 0, -10000, 10000);
RooRealVar b("b", "b", 0, -10000, 10000);
RooGenericPdf StraightLine("StraightLine", "Linear Model", "m*x + b", RooArgSet(x, m, b));
//StraightLine.fitTo(CombinedData,RooFit::Minimizer("Minuit","Minimize"),RooFit::Save(),NormSet(NormalizationSet));
StraightLine.fitTo(CombinedData,RooFit::Minimizer("Minuit","Minimize"),RooFit::Strategy(2),RooFit::Save(),NormSet(NormalizationSet));
// Extracting the parameters of the straight line fit
double P0 = m.getVal();
double P1 = b.getVal();
RooRealVar Slope("Slope", "Slope", P0, P0, P0);
RooRealVar YIntercept("YIntercept", "YIntercept", P1, P1, P1);
RooPolynomial Background("PolBackground", "First-Order Polynomial Background P.D.F.", Mass, RooArgList(Slope, YIntercept), 1);
// Fixing the parameters of the first-order polynomial background P.D.F.
Slope.setConstant(true);
YIntercept.setConstant(true);
RooRealVar SigAmplitude("SigAmplitude", "SigAmplitude", 50000, 0, 1000000);
RooRealVar BgAmplitude("BgAmplitude", "BgAmplitude", 10000, 0, 1000000);
RooExtendPdf ExtSignal("ExtSignal", "ExtSignal", Signal, SigAmplitude);
RooExtendPdf ExtBackground("ExtBackground", "ExtBackground", Background, BgAmplitude);
RooAddPdf* ModelPDF = new RooAddPdf("ModelPDF","ModelPDF",RooArgList(ExtSignal,ExtBackground));
ModelPDF->fixCoefNormalization(NormalizationSet);
//RooFitResult* FitResult = ModelPDF->fitTo(DataHist,NumCPU(4,kTRUE),Save(),Extended(),RooFit::Minimizer("Minuit","Minimize"),NormSet(NormalizationSet));
RooFitResult* FitResult = ModelPDF->fitTo(DataHist,NumCPU(4,kTRUE),Save(),Extended(),RooFit::Minimizer("Minuit","Minimize"),RooFit::Strategy(2),NormSet(NormalizationSet));
double SignalYield = SigAmplitude.getVal();
double SignalYieldError = SigAmplitude.errorVar()->getVal();
double BackgroundYield = BgAmplitude.getVal();
double BackgroundYieldError = BgAmplitude.errorVar()->getVal();
double TotalYield = SignalYield + BackgroundYield;
double NumberOfEvents = Data->sumEntries();
// Plotting Data onto Frame
TString FrameTitle;
TString XAxisTitle;
TString YAxisTitle;
FrameTitle = "#font[12]{B Mass Distribution}";
YAxisTitle = "#font[12]{Events / 2.5 MeV}";
XAxisTitle = "#font[12]{B^{#pm}#rightarrowK^{#pm}J/#it{#psi} Invariant Mass[GeV]}";
// Plotting Data Points
RooPlot* Frame = Mass.frame(Title(FrameTitle),Bins(NBins),Range(LeftCut,RightCut));
Data->plotOn(Frame,Name("Data"),XErrorSize(0));
double BgNormFactor = BgAmplitude.getVal()/SigAmplitude.getVal();
// Plotting Full Model
ModelPDF->plotOn(Frame,LineColor(kRed),LineWidth(3),NumCPU(1,kTRUE),Name("Model"));
ModelPDF->plotOn(Frame,Components("ExtBackground"),LineWidth(3),LineColor(kAzure),LineStyle(9),Name("Background"));
ModelPDF->plotOn(Frame,Components("ExtSignal"),LineWidth(3),DrawOption("F"),FillColor(kGreen),LineColor(kGreen),Name("Double-Gaussian"));
ModelPDF->plotOn(Frame,Components("ExtSignal"),LineWidth(3),DrawOption("L"),LineColor(kOrange),Name("GaussianBackground"),Normalization(BgNormFactor));
ModelPDF->plotOn(Frame,Components("Gaussian1"),LineWidth(3),LineColor(kSlate),Name("Gaussian1"));
ModelPDF->plotOn(Frame,Components("Gaussian2"),LineWidth(3),LineColor(kPurple),Name("Gaussian2"));
/// Creating and Filling Canvas with Plots and Info
TCanvas* MassProjection = new TCanvas("Mass Fit", "Mass Fit", 900, 900);
MassProjection->SetWindowSize(1920,1080);
MassProjection->cd();
Frame->GetXaxis()->SetTitleOffset(1.);
Frame->GetXaxis()->SetTitle(XAxisTitle);
Frame->GetYaxis()->SetTitle(YAxisTitle);
Frame->SetMinimum(1e-5);
Frame->Draw();
double Offset1 = 0.020;
double Offset2 = 0.161;
TString DecayMode;
if (SignalType == "J" || JorP == "J") {
if (FinalStateLeptons == "MuMu")
DecayMode = "#font[12]{B^{#pm}#rightarrowK^{#pm}J/#it{#psi}#rightarrow#it{#mu^{+}}#it{#mu^{-}}}";
else if (FinalStateLeptons == "EMu")
DecayMode = "#font[12]{B^{#pm}#rightarrowK^{#pm}J/#it{#psi}#rightarrowe^{#pm}#it{#mu^{#mp}}}";
}
else if (SignalType == "P" || JorP == "P") {
if (FinalStateLeptons == "MuMu")
DecayMode = "#font[12]{B^{#pm}#rightarrowK^{#pm}#it{#psi}(2S)#rightarrow#it{#mu^{+}}#it{#mu^{-}}}";
else if (FinalStateLeptons == "EMu")
DecayMode = "#font[12]{B^{#pm}#rightarrowK^{#pm}#it{#psi}(2S)#rightarrowe^{#pm}#it{#mu^{#mp}}}";
}
TLatex* Latex = new TLatex();
Latex->SetTextSize(0.03);
Latex->DrawLatex(LeftCut+Offset1, Frame->GetMaximum()*0.900, "#font[22]{#color[1]{#it{CMS Preliminary}}}");
Latex->DrawLatex(LeftCut+Offset1, Frame->GetMaximum()*0.850, "#font[22]{#it{#sqrt{s} = 13 TeV, L = 137.19 fb^{-1}}}");
Latex->DrawLatex(LeftCut+Offset1, Frame->GetMaximum()*0.800, DecayMode);
Latex->SetTextSize(0.02);
Latex->DrawLatex(LeftCut+Offset2, Frame->GetMaximum()*0.620,Form("N(Signal) = %.2f #pm %.2f", SignalYield, SignalYieldError));
Latex->DrawLatex(LeftCut+Offset2, Frame->GetMaximum()*0.570,Form("N(Background) = %.2f #pm %.2f", BackgroundYield, BackgroundYieldError));
Latex->DrawLatex(LeftCut+Offset2, Frame->GetMaximum()*0.520,Form("Total number of events returned from the fit = %.2f", TotalYield));
Latex->DrawLatex(LeftCut+Offset2, Frame->GetMaximum()*0.470,Form("Total number of events = %.0f", NumberOfEvents));
RooCurve* PlotData = Frame->getCurve("Data");
RooCurve* PlotModel = Frame->getCurve("Model");
RooCurve* PlotDG = Frame->getCurve("Double-Gaussian");
RooCurve* PlotBG = Frame->getCurve("Background");
RooCurve* PlotG1 = Frame->getCurve("Gaussian1");
RooCurve* PlotG2 = Frame->getCurve("Gaussian2");
TLegend* MassLeg = new TLegend(0.63, 0.69, 0.83, 0.85);
MassLeg->SetTextFont(62);
MassLeg->SetTextSize(0.02);
MassLeg->SetBorderSize(0);
MassLeg->SetFillColor(0);
MassLeg->AddEntry(PlotData, "Data", "lep");
MassLeg->AddEntry(PlotModel,"Total", "l");
MassLeg->AddEntry(PlotDG, "Signal", "f");
MassLeg->AddEntry(PlotBG, "Background", "l");
MassLeg->AddEntry(PlotG1, "First Gaussian Components", "l");
MassLeg->AddEntry(PlotG2, "Second Gaussian Components", "l");
MassLeg->Draw();
// Getting Value of chi2/ndf
double Chi2NDF = Frame->chiSquare("Model", "Data", FitResult->floatParsFinal().getSize());
double FitQuality = FitResult->covQual();
double FitStatus = FitResult->status();
FilePath = OutputDirectoryName + "/Fit.png";
MassProjection->SaveAs(FilePath.c_str());
FilePath.clear();
std::cout << endl;
std::cout << "#####################################################################" << endl;
std::cout << "Chi-Squared over Number of Degrees of Freedom = " << Chi2NDF << endl;
std::cout << endl;
std::cout << "Covariance (Fit) Quality = " << FitQuality << endl;
std::cout << endl;
std::cout << "Fit Status = " << FitStatus << endl;
std::cout << "#####################################################################" << endl;
std::cout << endl;
RooArgSet* Parameters = ModelPDF->getParameters(*Data);
TIterator* iPar = Parameters->createIterator();
RooRealVar* Parameter;
while ((Parameter = (RooRealVar*)iPar->Next())) {
std::cout << "Parameter: " << setw(20) << left << Parameter->GetName() << "\tValue: " << setw(10) << fixed << setprecision(6) << Parameter->getVal();
std::cout << "\tUncertainty: " << setw(10) << fixed << setprecision(6) << Parameter->getError();
TextFile << "Parameter: " << setw(20) << left << Parameter->GetName() << "\tValue: " << setw(10) << fixed << setprecision(6) << Parameter->getVal();
TextFile << "\tUncertainty: " << setw(10) << fixed << setprecision(6) << Parameter->getError();
if (Parameter->isConstant()) {
std::cout << " (Fixed)";
TextFile << " (Fixed)";
}
std::cout << endl;
TextFile << endl;
}
delete iPar;
delete Parameters;
std::cout << endl;
TextFile << endl;
// Significance of Signal
RooArgSet NormalizationSet(Mass);
Signal.fixCoefNormalization(NormalizationSet);
// Creating Likelihood Function
RooRealVar SigYield("SigYield", "Signal Yield", SigAmplitude.getVal());
RooRealVar BgYield("BgYield", "Background Yield", BgAmplitude.getVal());
RooFormulaVar NegLogLiklihood("NegLogLiklihood", "-log(SigYield * ExtSignal + BgYield * ExtBackground)", RooArgList(SigYield, BgYield, ExtSignal, ExtBackground));
// Creating Minuit2 minimizer and setting initial values of parameters
RooMinimizer Minimizer(NegLogLiklihood);
SigAmplitude.setVal(0.);
BgAmplitude.setVal(1.);
Minimizer.setMinimizerType("Minuit2");
Minimizer.minimize("Minuit2", "Migrad");
Minimizer.setPrintLevel(-1);
// Calculating Profile Likelihood Ratio
double NLLhMin = NegLogLiklihood.getVal(); // <---<< The value of the negative log-likelihood function evaluated at the minimum found by the minimizer
SigAmplitude.setVal(0.); // <---<< Fitting the model under the assumption that the signal component has zero yield
Minimizer.minimize("Minuit2", "Migrad");
double NLLhNull = NegLogLiklihood.getVal(); // <---<< The minimum value of the log-likelihood function under the null hypothesis
double ProfileLikelihoodRatio = 2*(NLLhNull-NLLhMin);
// Calculating Significance
double PValue = TMath::Prob(ProfileLikelihoodRatio, 94); // <---<< 94 degrees of freedom
double SignalSignificance = TMath::NormQuantile(1-PValue);
std::cout << "P-Value = " << PValue << std::endl;
std::cout << "Significance of Signal = " << SignalSignificance << std::endl;
TextFile << "P-Value = " << PValue << std::endl;
TextFile << "Significance of Signal = " << SignalSignificance << std::endl;
std::cout << endl;
TextFile << endl;
}
