Convolution fit issues,

Hirmans1 · April 7, 2025, 9:26am

Dear Expert,

I’m trying to fit the invariant mass with the square modulus of a complex function convoluted with a Gaussian function. However, when I add the Gaussian function to the complex function, the fit fails, as the graph shows. I’m not quite sure why the fit fails when using the complex function with a Gaussian.
Could you please provide me with a solution to this issue?

Thanks in advance.

Sorry for not providing the plots, cause I am new user. could you help me how to upload them.

here is the code im using :

include <TFile.h>

#include <TTree.h>

#include <TH1F.h>

#include <TF1.h>

#include <TCanvas.h>

#include <TLegend.h>

#include <TStyle.h>

#include <TGaxis.h>

#include <TLine.h>

#include <TMath.h>

#include <cmath>

#include <TF1Convolution.h>

#include <iostream>

#include <sstream>

#include <iomanip>

**using** **namespace** std;

// Function declarations

**double** Px(**double** x, **double** c0, **double** c1, **double** c2);

**double** pK(**double** x, **double** mK);

**double** pPhi(**double** x, **double** mJpsi, **double** mPi0);

**double** B(**double** r, **double** p);

**double** phs(**double** x, **double** mK, **double** M0);

**double** Gamma_m(**double** x, **double** M0, **double** Gamma0, **double** mK, **double** r);

**double** f_m(**double** x, **double** M0, **double** Gamma0, **double** mJpsi, **double** mK, **double** mPi0, **double** r);

**double** k(**double** M0, **double** Gamma0);

**double** bkg_nonRel(**double** x, **double** * par);

**double** res_nonRel(**double** x, **double** R, **double** M0, **double** Gamma0, **double** mK, **double** r);

**double** interference_nonRel(**double** x, **double** c0, **double** c1, **double** c2, **double** R, **double** delta, **double** M0, **double** Gamma0, **double** mK, **double** r);

**double** dN_dx(**double** x, **double** c0, **double** c1, **double** c2, **double** R, **double** delta, **double** M0, **double** Gamma0, **double** mK, **double** r);

**double** sqrtR_plus(**double** M0, **double** c0, **double** c1, **double** c2, **double** delta, **double** R, **double** Gamma0, **double** mK, **double** r);

**double** sqrtR_minus(**double** M0, **double** c0, **double** c1, **double** c2, **double** delta, **double** R, **double** Gamma0, **double** mK, **double** r);

**double** sig_res_rel(**double** m, **double** * par);

**double** Interference_res_rel(**double** m, **double*** par);

**double** FullPdf(**double** * x, **double** * par);

**double** Background(**double** * x, **double** * par);

**double** Signal(**double** * x, **double** * par);

**double** Interference(**double** * x, **double** * par);

**double** Gaussian(**double** * x, **double** * par);

// Function implementations

**double** Px(**double** x, **double** c0, **double** c1, **double** c2) {

**return** c0 + c1 * x + c2 * x * x;

}

**double** pK(**double** x, **double** mK) {

**double** val = x*x - 4.0*mK*mK;

**if** (val <= 0) **return** 0.0; // Return 0 instead of -1 for invalid momentum

**return** 0.5 * sqrt(val);

}

**double** pPhi(**double** x, **double** mJpsi, **double** mPi0) {

**double** term = (mJpsi*mJpsi - (x + mPi0)*(x + mPi0)) * (mJpsi*mJpsi - (x - mPi0)*(x - mPi0));

**if** (term <= 0) **return** 0.0;

**return** sqrt(term) / (2.0 * mJpsi);

}

**double** B(**double** r, **double** p) {

**if** (p <= 0) **return** 1.0; // Avoid division by zero

**return** 1.0 / sqrt(1.0 + (r * p) * (r * p));

}

**double** phs(**double** x, **double** mK, **double** M0) {

**double** pK_m = pK(x, mK);

**double** pK_M0 = pK(M0, mK);

**return** pK_m * pK_M0;

}

**double** f_m(**double** x, **double** M0, **double** Gamma0, **double** mJpsi, **double** mK, **double** mPi0, **double** r) {

**double** pK_m = pK(x, mK);

**double** pK_M0 = pK(M0, mK);

**double** pPhi_m = pPhi(x, mJpsi, mPi0);

**double** pPhi_M0 = pPhi(M0, mJpsi, mPi0);

**if** (pK_M0 == 0 || pPhi_M0 == 0) **return** 0.0;

**return** (pPhi_m / pPhi_M0) * (pK_m / pK_M0) * (B(r, pPhi_m) / B(r, pPhi_M0)) * (B(r, pK_m) / B(r, pK_M0));

}

**double** Gamma_m(**double** x, **double** M0, **double** Gamma0, **double** mK, **double** r) {

**double** pK_m = pK(x, mK);

**double** pK_M0 = pK(M0, mK);

**if** (pK_M0 == 0 || x == 0) **return** 0.0;

**return** pow(pK_m / pK_M0, 3) * (M0 / x) * Gamma0 * (B(r, pK_m) / B(r, pK_M0));

}

**double** k(**double** M0, **double** Gamma0) {

**double** gam = M0 * sqrt(M0*M0 + Gamma0*Gamma0);

**return** (2.0 * sqrt(2.0) * M0 * Gamma0 * gam) / (M_PI * (M0*M0 + gam));

}

// I-Non-resonance contribution:

**double** bkg_nonRel(**double** x, **double** * par) {

**return** Px(x, par[0], par[1], par[2]);

}

**double** res_nonRel(**double** x, **double** R, **double** M0, **double** Gamma0, **double** mK, **double** r) {

**double** gamma_m = Gamma_m(x, M0, Gamma0, mK, r);

**if** (gamma_m <= 0) **return** 0.0;

**return** R * gamma_m / (2.0 * M_PI) / ((x - M0)*(x - M0) + gamma_m*gamma_m/4.0);

}

**double** interference_nonRel(**double** x, **double** c0, **double** c1, **double** c2, **double** R, **double** delta, **double** M0, **double** Gamma0, **double** mK, **double** r) {

**double** P_m = Px(x, c0, c1, c2);

**double** res = res_nonRel(x, R, M0, Gamma0, mK, r);

**return** P_m + res * cos(delta);

}

**double** dN_dx(**double** x, **double** c0, **double** c1, **double** c2, **double** R, **double** delta, **double** M0, **double** Gamma0, **double** mK, **double** r) {

**double** params[3] = {c0, c1, c2}; // Avoid memory leak

**double** bkg = bkg_nonRel(x, params);

**double** res = res_nonRel(x, R, M0, Gamma0, mK, r);

**double** interfer = interference_nonRel(x, c0, c1, c2, R, delta, M0, Gamma0, mK, r);

**return** fmax(bkg + res + interfer, 0.0); // Prevent negative values

}

// Two solutions from x = M0

**double** sqrtR_plus(**double** M0, **double** c0, **double** c1, **double** c2, **double** delta, **double** R, **double** Gamma0, **double** mK, **double** r) {

**double** params[3] = {c0, c1, c2};

**double** gamma_M0 = Gamma_m(M0, M0, Gamma0, mK, r);

**double** P_M0 = bkg_nonRel(M0, params);

**double** dN_dx_M0 = dN_dx(M0, c0, c1, c2, R, delta, M0, Gamma0, mK, r);

**double** discriminant = (M_PI * gamma_M0 * P_M0 / 2.0) * pow(sin(delta), 2) - (M_PI * gamma_M0 / 2.0) * (P_M0 - dN_dx_M0);

**if** (discriminant < 0) **return** 0.0; // Handle negative discriminant

**return** -sqrt(M_PI * gamma_M0 * P_M0 / 2.0) * sin(delta) + sqrt(discriminant);

}

**double** sqrtR_minus(**double** M0, **double** c0, **double** c1, **double** c2, **double** delta, **double** R, **double** Gamma0, **double** mK, **double** r) {

**double** params[3] = {c0, c1, c2};

**double** gamma_M0 = Gamma_m(M0, M0, Gamma0, mK, r);

**double** P_M0 = bkg_nonRel(M0, params);

**double** dN_dx_M0 = dN_dx(M0, c0, c1, c2, R, delta, M0, Gamma0, mK, r);

**double** discriminant = (M_PI * gamma_M0 * P_M0 / 2.0) * pow(sin(delta), 2) - (M_PI * gamma_M0 / 2.0) * (P_M0 - dN_dx_M0);

**if** (discriminant < 0) **return** 0.0;

**return** -sqrt(M_PI * gamma_M0 * P_M0 / 2.0) * sin(delta) - sqrt(discriminant);

}

// II-Resonance contribution:

**double** sig_res_rel(**double** m, **double** * par) {

**double** c0 = par[0];

**double** c1 = par[1];

**double** c2 = par[2];

**double** M0 = par[3];

**double** Gamma0 = par[6], mJpsi = par[7], mK = par[8];

**double** mPi0 = par[9], r = par[10];

**double** delta = par[5];

**double** R = par[4];

**double** k_par = k(M0, Gamma0);

**double** sqrtK = sqrt(k_par);

**double** gamma_m = Gamma_m(m, M0, Gamma0, mK, r);

**double** f = f_m(m, M0, Gamma0, mJpsi, mK, mPi0, r);

**double** ReR = sqrtR_plus(M0, c0, c1, c2, delta, R, Gamma0, mK, r) * cos(delta);

**double** ImR = sqrtR_plus(M0, c0, c1, c2, delta, R, Gamma0, mK, r) * sin(delta);

**double** Aphi_squared_mod_m = (ReR * ReR + ImR * ImR) / ((m * m - M0 * M0) * (m * m - M0 * M0) + m * m * gamma_m * gamma_m);

**return** sqrtK * sqrtK * f * f * Aphi_squared_mod_m * phs(m, mK, M0);

}

**double** Interference_res_rel(**double** m, **double*** par) {

**double** bkg = bkg_nonRel(m, par);

**double** phs_x = phs(m, par[8], par[3]);

**double** gamma_m = Gamma_m(m, par[3], par[6], par[8], par[10]);

**double** k_par = k(par[3], par[6]);

**double** sqrtK = sqrt(k_par);

**double** f = f_m(m, par[3], par[6], par[7], par[8], par[9], par[10]);

**double** ReR = sqrtR_plus(par[3], par[0], par[1], par[2], par[5], par[4], par[6], par[8], par[10]) * cos(par[5]);

**double** ImR = sqrtR_plus(par[3], par[0], par[1], par[2], par[5], par[4], par[6], par[8], par[10]) * sin(par[5]);

**return** 2.0 * sqrt(bkg * phs_x) * sqrtK * f * (ReR * (m * m - par[3] * par[3]) + ImR * (m * gamma_m)) /

((m * m - par[3] * par[3]) * (m * m - par[3] * par[3]) + m * m * gamma_m * gamma_m);

}

// Model

**double** FullPdf(**double** * x, **double** * par) {

**double** m = x[0];

**double** bkg = bkg_nonRel(m, par);

**double** sig_res = sig_res_rel(m, par);

**double** Inter = Interference_res_rel(m, par);

**return** bkg + sig_res + Inter;

}

**double** Background(**double** * x, **double** * par) {

**return** bkg_nonRel(x[0], par);

}

**double** Signal(**double** * x, **double** * par) {

**return** sig_res_rel(x[0], par);

}

**double** Interference(**double** * x, **double** * par) {

**return** Interference_res_rel(x[0], par);

}

// Gaussian function:

**double** Gaussian(**double** * x, **double** * par) {

**double** m = x[0];

**double** mu = par[0];

**double** sigma_m = par[1];

**return** 1.0 / (sigma_m * sqrt(2 * M_PI)) * exp(- pow((m - mu), 2) / (2 * pow(sigma_m, 2)));

}

/*double AddPdf(double *x, double * par){

double m = x[0];

double f1 = FullPdf(&m, par);

double f2 = Gaussian(&m, par);

return f1 * f2;

}

*/

**int** main() {

gStyle->SetOptFit(1111);

// Open the ROOT file

**const** **char*** homeDir = getenv("HOME");

TString filePath = TString::Format("%s/Desktop/Jpsi_211_cut1.root", homeDir);

TFile* file = TFile::Open(filePath);

**if** (!file || file->IsZombie()) {

cerr << "Error opening the file!" << endl;

**return** 1;

}

// Get the TTree

TTree* tree = (TTree*)file->Get("Jpsi_211");

**if** (!tree) {

cerr << "TTree not found!" << endl;

**return** 1;

}

// Variable to hold mphi

Double_t mphi;

cout << "mphi: " << mphi << endl;

**if** (tree->SetBranchAddress("mphi", &mphi) != 0) {

cerr << "Error: Could not set branch address!" << endl;

**return** 1;

}

// Create a histogram for mphi

TH1F* hist = **new** TH1F("fit_result", "Background + Signal", 100, 0.98, 1.18);

// Fill the histogram

Long64_t nEntries = tree->GetEntries();

**for** (Long64_t i = 0; i < nEntries; i++) {

tree->GetEntry(i);

hist->Fill(mphi);

}

// Check if no entries were found

**if** (nEntries == 0) {

cerr << "Warning: No entries were found in the TTree!" << endl;

}

// Create a canvas and draw the histogram

TCanvas* canvas = **new** TCanvas("canvas", "Canvas", 800, 600);

hist->SetTitle("3097MeV");

hist->GetXaxis()->SetTitle("M(K^{+}K^{-}) [GeV/c^{2}]");

hist->GetXaxis()->SetNdivisions(510);

hist->GetXaxis()->SetLabelSize(0.04);

// Set the Y-axis title

ostringstream yAxisTitle;

yAxisTitle << "Events / " << fixed << setprecision(2) << (hist->GetBinWidth(1) * 1000) << " MeV/c^{2}";

hist->GetYaxis()->SetTitle(yAxisTitle.str().c_str());

hist->GetYaxis()->SetRangeUser(-1400, hist->GetMaximum() * 1.2);

hist->SetMarkerStyle(kFullCircle);

hist->Draw("E1 PMC");

// Draw the horizontal line at y = 0

TLine* line = **new** TLine(hist->GetXaxis()->GetXmin(), 0, hist->GetXaxis()->GetXmax(), 0);

line->SetLineColor(kBlack);

line->Draw("SAME");

// Define the full PDF

TF1* fullPDF = **new** TF1("fullPDF", FullPdf, 0.98, 1.18, 11);

fullPDF->SetParameters(0.001, -0.0001, -0.002, 1.019461, 3.1, -M_PI/4.0, 0.00426, 3.0969, 0.493677, 0.1349768, 3.0);

// Set limits for parameters

fullPDF->SetParLimits(0, -100000, 100000); // c0

fullPDF->SetParLimits(1, -100000, 100000); // c1

fullPDF->SetParLimits(2, -80000, 80000); // c2

fullPDF->SetParLimits(3, 1.019461 , 1.019461); // phi mass limits (M0)

fullPDF->SetParLimits(4, .0001, 1000); // R limits

fullPDF->SetParLimits(5, -M_PI, 0); // Delta limits

fullPDF->SetParLimits(6, 0.004249, 0.004249); // Gamma0 limits

fullPDF->SetParLimits(7, 3.0969, 3.0969); // J/psi mass limits (mJpsi)

fullPDF->SetParLimits(8, 0.493677, 0.493677); // Kaon mass limits (mK)

fullPDF->SetParLimits(9, 0.1349768,0.1349768 ); // Neutral pion mass (mPi0) (0.1349768)

fullPDF->SetParLimits(10, 3.0, 3.0); // r (fixed)

fullPDF->SetParNames("c0", "c1", "c2", "M_{#phi}", "R_{1}", "#delta_{1}", "#Gamma_{0}", "M_{J/#psi}", "M_{K}", "M_{#pi^{0}}", "r");

// Fit the histogram with the full PDF

**int** fitStatus = hist->Fit(fullPDF, "R");

**if** (fitStatus != 0) {

cerr << "Fit failed with status: " << fitStatus << endl;

}**else** {cerr <<"Fit success with status:"<< fitStatus <<endl;}

// Draw the fitted function

fullPDF->SetLineColor(kBlue-3);

fullPDF->Draw("SAME");

// Define the background

TF1 *background = **new** TF1("background", Background, 0.98, 1.18, 3);

background->SetParameters(fullPDF->GetParameters());

background->SetLineColor(kBlue);

background->SetLineStyle(2);

background->Draw("SAME");

TF1* signal = **new** TF1("signal", Signal, 0.98, 1.18, 11);

signal->SetParameters(fullPDF->GetParameters());

signal->SetLineColor(kRed);

signal->SetLineStyle(3);

signal->Draw("SAME");

TF1* interference = **new** TF1("interference", Interference, 0.98, 1.18, 11);

interference->SetParameters(fullPDF->GetParameters());

interference->SetLineColor(kGreen);

interference->SetLineStyle(4);

interference->Draw("SAME");

TF1 *gaussian = **new** TF1("gaussian", Gaussian, 0.98, 1.18, 2);

gaussian->SetParameter(0, 1.0); // Mean

gaussian->SetParameter(1, 0.01); // Sigma (example value)

// Create the convolution of the full PDF and the Gaussian

TF1Convolution* conv = **new** TF1Convolution(fullPDF, gaussian, 0.98, 1.18, **true**);

conv->SetRange(0.98, 1.18);

TF1* convFunc = **new** TF1("convFunc", conv, 0.98, 1.18, conv->GetNpar());

**double** convParams[13]; // 11 from fullPDF + 2 from Gaussian

**for** (**int** i = 0; i < 11; i++)

convParams[i] = fullPDF->GetParameter(i);

convParams[11] = gaussian->GetParameter(0); // Gaussian mean

convParams[12] = gaussian->GetParameter(1); // Gaussian sigma

convFunc->SetParameters(convParams);

convFunc->SetParameters(fullPDF->GetParameters());

//hist ->Fit(convFunc, "R");

// convFunc -> SetLineColor(kBlue);

//convFunc->Draw("same");

// Create a legend to show the models

TLegend* legend = **new** TLegend(0.15, 0.738, 0.3, 0.89);

legend->AddEntry(hist, "Data", "ep");

legend->SetBorderSize(0);

legend->AddEntry(fullPDF, "Total fit", "l");

legend->AddEntry(background, "Comp1 (Bkg)", "l");

legend->AddEntry(signal, "Comp2 (Sig)", "l");

legend->AddEntry(interference, "Comp3 (Interference)", "l");

legend->Draw();

// Save the canvas as a PDF

canvas->SaveAs("~/Desktop/Plots/fit_dat0919_p_c2.pdf");

// Clean up

**delete** canvas;

**delete** hist;

**delete** fullPDF;

**delete** gaussian;

**delete** convFunc;

**delete** conv;

**delete** background;

**delete** signal;

**delete** interference;

**return** 0;

}

devajith · April 7, 2025, 9:31am

Hi @Hirmans1,

Thank you for your question and welcome to the ROOT forum.

Can you give us a minimal reproducer? That will help us pinpoint your issue better.
In the meantime, I’m adding our RooFit expert @jonas in the loop.

Cheers,
Dev

Hirmans1 · April 8, 2025, 7:40am

Okay, but I have to upload some info, cause the code that I pasted here yesterday changed. May be I have to upload it.

Thanks,
Hirmans/