# Fit using Double sided crystal ball function

Dear root experts,

I implemented a double-sided crystal ball function to be used to fit a histogram.
I based on this article page 10 to implement my code: https://dash.harvard.edu/bitstream/handle/1/29362185/1606.03833.pdf?sequence=1
The code seems to work but I am not satisfied with the fit results, I was expected to be able to fit the right tail of the histogram for instance given the high number of parameters that one can tune to adjust the fit. So I was wondering if someone more experienced on fitting could help.
I am attaching the code, the input histogram and the result I got.
Thanks.
The code can be run by just: root -l tmp.C

DSCBFit.pdf (20.4 KB) signal_50_400_Nominal_reco_m4l_output_combined.root (382.7 KB) ___

``````double DoubleSidedCrystalballFunction(double *x, double *par)
{
double alpha_l = par;
double alpha_h = par;
double n_l     = par;
double n_h     = par;
double mean	= par;
double sigma	= par;
double N	= par;
float t = (x-mean)/sigma;
double result;
double fact1TLessMinosAlphaL = alpha_l/n_l;
double fact2TLessMinosAlphaL = (n_l/alpha_l) - alpha_l -t;
double fact1THihgerAlphaH = alpha_h/n_h;
double fact2THigherAlphaH = (n_h/alpha_h) - alpha_h +t;

if (-alpha_l <= t && alpha_h >= t)
{
result = exp(-0.5*t*t);
}
else if (t < -alpha_l)
{

result = exp(-0.5*alpha_l*alpha_l)*pow(fact1TLessMinosAlphaL*fact2TLessMinosAlphaL, -n_l);

}
else if (t > alpha_h)
{
result = exp(-0.5*alpha_l*alpha_l)*pow(fact1THihgerAlphaH*fact1THihgerAlphaH, -n_h);

}
return N*result;
}

void tmp()
{
TCanvas *c1 = new TCanvas("c1","c1",700,500);
TFile *f = new TFile("signal_50_400_Nominal_reco_m4l_output_combined.root");
TH1F *hpx = (TH1F*)f->Get("havgM_all");
TF1 *fitDSCB = new TF1("fitDSCB",DoubleSidedCrystalballFunction, 255,420, 7);
hpx->GetYaxis()->SetRangeUser(0, 30);
gStyle->SetOptFit(1);
fitDSCB->SetParameters(1, 2, 2, 1, hpx->GetMean(), hpx->GetRMS(),hpx->Integral(255, 425));
fitDSCB->SetParNames ("alpha_{low}","alpha_{high}","n_{low}", "n_{high}", "mean", "sigma", "Norm");
hpx->Fit(fitDSCB, "", "", 255, 420);
Double_t chi2 = fitDSCB->GetChisquare();
cout << "chi2 = " << chi2 << endl;
hpx->Draw("func");
hpx->GetXaxis()->SetRangeUser(200, 550);

}

``````

ROOT Version: 6.18/04

Built for macosx64 on Sep 11 2019, 15:38:23

From tags/v6-18-04@v6-18-04

1 Like

1 Like

It seems your fit function is not able to describe your data … it seems to me that you should modify at least one line:
`result = exp(-0.5*alpha_h*alpha_h)*pow(fact1THihgerAlphaH*fact2THigherAlphaH, -n_h);`

1 Like

Hi:

Yes, well spotted. I will fix that and try to fit it again.
That was the issue actually. The fit is better now.
I attached it. DSCBFit.pdf (20.4 KB)

Thanks.

Best Regards,

Diallo.

You could also try the `"L"` fit option.

1 Like

Hi:

I tried it, I can go even in a higher range.
But I also changed the parameters, so I am not sure if it’s only related to the “L” option or/and also the parameters I changed. DSCBFit.pdf (20.4 KB)

Thanks.

Best Regards,

Diallo

Hello,

I am not sure if I should open a new thread or not, but since it’s the same code I think it should be fine. Actually, if I tried to generate random histogram using the function I implemented, the code works if I use this for instance:

``````double DoubleSidedCrystalballFunction(double *x, double *par)
{
double alpha_l = par;
double alpha_h = par;
double n_l     = par;
double n_h     = par;
double mean	= par;
double sigma	= par;
double N	= par;
float t = (x-mean)/sigma;
double result;
double fact1TLessMinosAlphaL = alpha_l/n_l;
double fact2TLessMinosAlphaL = (n_l/alpha_l) - alpha_l -t;
double fact1THihgerAlphaH = alpha_h/n_h;
double fact2THigherAlphaH = (n_h/alpha_h) - alpha_h +t;

if (-alpha_l <= t && alpha_h >= t)
{
result = exp(-0.5*t*t);
}
else if (t < -alpha_l)
{

result = exp(-0.5*alpha_l*alpha_l)*pow(fact1TLessMinosAlphaL*fact2TLessMinosAlphaL, -n_l);

}
else if (t > alpha_h)
{
result = exp(-0.5*alpha_h*alpha_h)*pow(fact1THihgerAlphaH*fact2THigherAlphaH, -n_h);

}
return N*result;
}

void ForumRoot()
{

double M_S_j;
TF1 *rand_DSCB_mS = new TF1("rand_DSCB_mS",DoubleSidedCrystalballFunction, 120, 1500, 7);
TH1D* hGen_mS = new TH1D("hGen_mS" ,";m_{S} GeV;events/[2.5 GeV]", 552,120,1500);
for (int j = 0; j < 58; j++)
//for (int j = 0; j < 281; j++)

{

rand_DSCB_mS->SetParameters(1, 2, 1, 2, 249, 4, 1);
M_S_j = rand_DSCB_mS->GetRandom();
hGen_mS->Fill(M_S_j);
}
// hGen_mS->GetXaxis()->SetRangeUser(0, 300);

//hGen_mS->Draw("l");

}
``````

If I change this:

``````rand_DSCB_mS->SetParameters(1, 2, 1, 2, 249, 4, 1);
``````

To this:

``````rand_DSCB_mS->SetParameters(1, 2, 1, 2, 249, 5, 1);
``````

where you can see only ‘4’ has been replaced by ‘5’ I start having so many errors like that:

``````Error in <GSLError>: Error 18 in qags.c at 548 : cannot reach tolerance because of roundoff error
Warning in <TF1::IntegralOneDim>: Error found in integrating function rand_DSCB_mS in [216.600000,223.500000] using AdaptiveSingular. Result = 0.726260 +/- 0.000000  - status = 18
Info in <TF1::IntegralOneDim>: 		Function Parameters = { p0 =  1.000000 , p1 =  2.000000 , p2 =  1.000000 , p3 =  2.000000 , p4 =  249.000000 , p5 =  5.000000 , p6 =  1.000000 }
``````

Does someone know please why I am having such kind of errors and how to fix them?
Thanks

At least: `double t`

BTW. In ROOT, before you try to fit your graph or histogram, you MUST set “reasonable” initial values for ALL parameters of your function (except for some “built-in” formulas, like “gaus”, for which the standard fit procedure can automatically “guess” them), otherwise the fitting procedure may easily misbehave.

1 Like

Yes, I agree with the fit, the parameters have to be set reasonable. Actually the fit is not a problem anymore. The last issue I was having is those errors I showed when trying to generate randomly a histogram. But apparently just changing:

``````float t
``````

to

``````double t
``````

fixed it. Thanks very much.

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.