Fit using Double sided crystal ball function

ROOT
1

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[0]; 
   double alpha_h = par[1]; 
   double n_l     = par[2]; 
   double n_h     = par[3]; 
   double mean	= par[4]; 
   double sigma	= par[5];
   double N	= par[6];
   float t = (x[0]-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);
  TPad* thePad = (TPad*)c1->cd();
  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);
  thePad->Print("DSCBFit.pdf");
  
}

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
2

I think @moneta can help you.

1 Like
3

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
4

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.

5

You could also try the "L" fit option.

1 Like
6

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

7

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[0]; 
  double alpha_h = par[1]; 
  double n_l     = par[2]; 
  double n_h     = par[3]; 
  double mean	= par[4]; 
  double sigma	= par[5];
  double N	= par[6];
  float t = (x[0]-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

8

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
9

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.

10

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