RooCurve::average(Model) invalid range

Dear experts,

I have a Roo fitting code that works pretty good. By that, I mean the fit looks good and the yields are reasonable. However, in the printsouts of the code I see the following errors:

ERROR:InputArguments – RooCurve::average(Model) invalid range (5.16125,5.16125)
ERROR:InputArguments – RooCurve::average(Model) invalid range (5.16375,5.16375)
[#0] ERROR:InputArguments – RooCurve::average(Model) invalid range (5.16625,5.16625)
[#0] ERROR:InputArguments – RooCurve::average(Model) invalid range (5.16875,5.16875)

ERROR:InputArguments – RooCurve::average(Model) invalid range (5.39875,5.39875)

I don’t use the “RooCurve::average” explicitly in my code.

Another problem I have is the following warning, despite using “RooFit::Minimizer(“Minuit”);” in my code:
Warning in ROOT::Math::FitConfig::CreateMinimizer: Could not create the Minuit2 minimizer. Try using the minimizer Minuit
----> trying with improve

Are these errors important and how I can fix them?
Thank you in advance


Please read tips for efficient and successful posting and posting code

ROOT Version: Not Provided
Platform: Not Provided
Compiler: Not Provided


Dear @Natilus ,

Thanks for reaching out to the forum! Let me ask @jonas to guide you with your questions.

Cheers,
Vincenzo

Dear @vpadulan ,
It seems like @jonas is not available. Is there anyone else who can help me with my problem?
Thank you,

Hello @Natilus, thanks for the ping! If you have questions about RooFit, please open them in the “RooFit” category so I can see them better. I changed the category of this thread now so it stays on my radar.

Can you please fill out the information in the bottom of your post, about ROOT version et cetera? The minimizer problem is most likely depending on how you get ROOT.

And for the error with the range: all the ranges are empty (the interval is zero). I guess that makes them invalid? Can you maybe share your code so I can understand where these calls to RooCurve::average() come from?

Cheers,
Jonas

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;
}

Dear @jonas ,
Did you have any chance to look at my code?
Thanks

Hi @Natilus!

Thanks for sharing parts of your code. I can’t run it though, which is a problem for fully understanding it, as I can’t experiment with it. The header and main function is missing. It would have helped also to add the terminal output and resulting plot.

But just from reading it: I have a few suggestions to address your issues.

  1. About the large chi2: it seems that you’re using RooCurve::chiSquare(). I would not recommend this, for reasons I have explained in this post:
    How to correctly extract the chi2, ndf, p-value of a RooFitResult - #3 by astauffe
    Can you try with RooAbsPdf::createChi2()?

  2. It’s true, when you call RooAddPdf::getVal(normSet) you’ll get the normalized value, which is not matching with your unnormalized histograms. You can use RooAddPdf::expectedEvents() to get the post-fit value of the total yield, which would be this scaling factor.

Does this help? Sorry for the long latency in replies. We’re wrapping up the next ROOT release, so the forum is falling a bit behind :frowning:

Cheers,
Jonas

Hi @jonas ,
Thank you very much for your help.

I am happy to share my entire code including the header file. Can I just share it here?
I will try your suggestions and let you know soon.

Please let me know if/how I can share my code?

Thanks again for your help!