Piecewise function fits of TH1F

ROOT
1

Hi all,

I am trying to fit a combination of TF1’s to a ratio represented by a TH1F to make a piecewise fit parameterising this ratio for later use (an example of the ratio and fits in attachment).
fit.pdf (29.6 KB)

I tried using three separate TF1’s defined over three separate intervals to piece-wise define the various parts of the TH1F, but this naturally has discontinuities at the boundaries.

I thought the way around this would be to define another function which was the combination of all three following $ROOTSYS/tutorials/fit/multifit.C but this fit goes completely crazy (black curve in attachment) despite the individual ones being “reasonable” so obviously something went wrong, but I can’t figure out what the problem is:

   TF1* ratiofit_lo = new TF1("ratiofit_lo","pol3",
                                           ratio->GetBinLowEdge(2),
                                           ratio->GetBinLowEdge(ratio->FindBin(45.)));
    TF1* ratiofit_mid = new TF1("ratiofit_mid","pol3",
                           ratio->GetBinLowEdge(ratio->FindBin(45.)),
                           ratio->GetBinLowEdge(ratio->FindBin(200.)));
    TF1* ratiofit_hi = new TF1("ratiofit_hi","pol0",
                          ratio->GetBinLowEdge(ratio->FindBin(200.)),
                          ratio->GetBinLowEdge(ratio->GetNbinsX()+1));
   
    ratiofit_total = new TF1("ratiofit_total","ratiofit_lo+ratiofit_mid+ratiofit_hi",
                        ratio->GetBinLowEdge(2),ratio->GetBinLowEdge(ratio->GetNbinsX()+1));

    ratio->Draw("BE");
    ratio->Fit("ratiofit_lo","FEBR");
    ratio->Fit("ratiofit_mid","FEBR+");
    ratio->Fit("ratiofit_hi","FEBR+");
    
    Double_t tf_par[9];
    ratiofit_lo->GetParameters(&tf_par[0]);
    ratiofit_mid->GetParameters(&tf_par[4]);
    ratiofit_hi->GetParameters(&tf_par[8]);

    ratiofit_total->SetParameters(tf_par);
    
    ratio->Fit("ratiofit_total","FER+");

    ratiofit_lo->Draw("same");
    ratiofit_mid->Draw("same");
    ratiofit_hi->Draw("same");
    ratiofit_total->Draw("same");

Example output of the 1) lo, 2) mid, 3) hi and 4) total fits:

 FCN=0.129737 FROM MINOS     STATUS=SUCCESSFUL    112 CALLS         278 TOTAL
                     EDM=5.10156e-15    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                  PARABOLIC         MINOS ERRORS        
  NO.   NAME      VALUE            ERROR      NEGATIVE      POSITIVE   
   1  p0           6.81594e+00   2.88358e-02  -2.88197e-02   2.88197e-02
   2  p1          -1.65216e-02   5.43298e-03  -5.43065e-03   5.43065e-03
   3  p2          -7.94212e-04   2.93527e-04  -2.93454e-04   2.93454e-04
   4  p3           1.56027e-05   4.70869e-06  -4.70829e-06   4.70829e-06

 FCN=1.30009 FROM MINOS     STATUS=SUCCESSFUL    112 CALLS         270 TOTAL
                     EDM=1.49407e-14    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                  PARABOLIC         MINOS ERRORS        
  NO.   NAME      VALUE            ERROR      NEGATIVE      POSITIVE   
   1  p0           5.97036e+00   1.65204e-01  -1.65426e-01   1.65426e-01
   2  p1          -4.49985e-03   6.34583e-03  -6.35455e-03   6.35455e-03
   3  p2           3.14528e-05   7.46964e-05  -7.47989e-05   7.47989e-05
   4  p3          -1.26909e-07   2.70417e-07  -2.70778e-07   2.70778e-07

 FCN=0.0328523 FROM MINOS     STATUS=SUCCESSFUL      2 CALLS          22 TOTAL
                     EDM=1.25136e-24    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                  PARABOLIC         MINOS ERRORS        
  NO.   NAME      VALUE            ERROR      NEGATIVE      POSITIVE   
   1  p0           5.50545e+00   6.48449e-02  -6.48449e-02   6.48449e-02

 FCN=316.245 FROM MINOS     STATUS=PROBLEMS     4597 CALLS        5275 TOTAL
                     EDM=5.76097e-14    STRATEGY= 1  ERROR MATRIX UNCERTAINTY 100.0 per cent
  EXT PARAMETER                  PARABOLIC         MINOS ERRORS        
  NO.   NAME      VALUE            ERROR      NEGATIVE      POSITIVE   
   1  p0           2.30241e+00   7.65860e-04  -1.92511e-03   1.92511e-03
   2  p1          -1.83096e-02   3.83123e-05  -2.14398e-04   2.14398e-04
   3  p2           1.61843e-04   4.26920e-07  -3.38016e-06   3.38016e-06
   4  p3          -4.14801e-07   2.55799e-09  -1.23457e-08   1.23457e-08
   5  p4           5.97036e+00   1.41421e+00                            
   6  p5          -4.49985e-03   1.41421e+00                            
   7  p6           3.14528e-05   1.41421e+00                            
   8  p7          -1.26909e-07   1.41421e+00                            
   9  p8           5.50545e+00   1.41421e+00

So in the total fit parameters p4 through p8 are nonsense?

Does anyone have any idea what is going wrong here? – why the red/blue/green curves (ratiofit_lo mid and hi) are ok, but what should (from my understanding) be a continuous fit across the entire range with parameters taken from the three individual fits goes wrong?
[ Or alternatively what a better way of parameterising this TH1F is? ]

Thanks for any advice!
Darren

2

Hi,

when you are combine a pol3 + pol3 + a pol0 this is equivalent to just fitting a simple pol3. You have therefore only 4 independent parameters, the others do not make sense.
If you want to fit globally the 3 individual pieces, taking into account their ranges is better creating the combined function using some C++ code (i.e. creating TF1 using root.cern.ch/root/htmldoc/TF1.html#F3 or D or E) and not bby combining the formulas.

Lorenzo

3

Hi Lorenzo,

Thanks for your help, I see now my misunderstanding was that I thought the x-ranges of the individual functions would be preserved and thus the total fit would not devolve to a simple polynomial fit across the whole range.

(For those who come to this thread looking for answers) I simply fixed this issue by instead defining my own piecewise function along the following lines:

Double_t fit_pt(Double_t *x, Double_t *par) { Double_t xv =x[0]; Double_t f = (xv<40)*(par[0]+xv*par[1]+xv*xv*par[2]+xv*xv*xv*par[3]) + (xv>=40 && xv<=100)*(par[4]+xv*par[5]+xv*xv*par[6]+xv*xv*xv*par[7]) + (xv>100)*(par[8]+xv*par[9]+xv*xv*par[10]+xv*xv*xv*par[11]); return f; }

Thanks again for your help,
Darren

2 Likes