Setting Parameter Limits, Fixing Parameters and Their effects on fitting

Hi everybody,

I want to fit a polN function on a histogram in a given range. Let’s say this function is the 1st order pol1.
I have checked the fit options that I can use in the link. I am not using predefined function in the meantime.

That’s what I picked for the fitting options: h1->Fit(new_Bcg, “LRS”); I am not sure to add C and F as an option. I use “ROOT::Math::MinimizerOptions::SetDefaultMinimizer(“Minuit2”);” indeed in the very beginning of the code.

If I FIX the parameters as the values of what they has to be precisely , then the fit is nice. However, I can’t get the integral error of the fit in a certain range due to the FIXING.

On the other hand, if I set par. limits or just initialize them as to the precise values, the fit is terrible because of FLEXıBLE parameter CHOOSING.

What do you suggest?

That’s the line to read out some values after fit to give you an idea:

cout << “\nBackground integral: " << new_Bcg->Integral( Bcg_x1_LowE , Bcg_x2_UpE) / BinWidth <<” ± "
<< new_Bcg->IntegralError( Bcg_x1_LowE , Bcg_x2_UpE, p_bcg, cov_mat.GetMatrixArray()) / BinWidth << endl ;


Without fixing the parameters. OBVIOUSLY THE BAD ONE.


After fixing the parameters. CORRECT ONE.

You cannot fit the “background” alone. You need to fit with a common “signal + background” function.

Normally, that’s what we do, but I am trying new stuff in my project. If I summarize my purpose by referencing the 2nd picture, I’ll draw a line in 2 coordinates by using a function , basically the pol1 at this moment. Then, I need 2 information from that function: 1st one is the integral and 2nd one is the error within some range I choose. The specific reason why I thought the function might be the good choice was that I would need to get an integral and its error in a smaller range on that drawn line.

I might be wrong. All these correlations and covaraiance matrix stuff to find the error from a TF1 object might need a total function. Is that what you are saying? I didn’t want to create a representating background histogram instead of this line looking like a graph, so I can use IntegralAndError() method ,instead.

Can you suggest and lead me the right direction by telling your idea?
Are you suggesting me to fit the total function as usual but just to use the background information in it?
I notice that if I use total function to fit and follow the normal steps, I get the different results when I compare to get those integrals from directly histogram bin content as TH1 object. Bin width is 1, so it’s no the issue. On the other hand, if I follow the same steps by fixing the parameters since I know what exactly they are, I get the SAME RESULTS as what I got from histogram.

Cheers.

Thanks for the old post. That’s my real sample data indeed.

Anyways. my code is already similar for these parts. However, you and me, like others, always focus on the net area and its error. The real problem is the estimation of a background indeed.

When I quickly ckeck the code, you emphasized the “fixed background” part in these lines:

if (fit_with_fixed_background) {
f1->FixParameter(3, f1->GetParameter(3));
f1->FixParameter(4, f1->GetParameter(4));
g1->Fit(f1, “WEMR”);
std::cout << “note: background parameters are fixed” << std::endl;
std::cout << "gaus peak integral = " << f1->GetParameter(0)
<< " +/- " << f1->GetParError(0)
<< std::endl;
}

Here again, you comment out the gauss integral while fixing the background parameters. Since the fixed parameters error will be ZERO, then it should also affect the result of gauss integral’s error.

I wonder what happens if you print out the background integral and its error in your code!

Does using graph instead of histogram or using gausn instead of gaus make ANY DIFFERENCE on results?

You can get the integral and error of the “signal + background” directly from the fitted “f1”. For the “signal” alone, you then have “fg”. If you want the “background” alone, then you will need to create another function (using “pol1”) and play the same game as for the “fg”.

BTW. Of course, if you fix the “background” parameters when fitting, then the “background” integral’s error will be zero because the corresponding covariance matrix elements will be zero. Obviously, the errors of the “signal” parameters will be different, too.

I know what you mean, we have done this many times. I assume you are referring this part of your code I show below:


TF1 *fg = new TF1(“fg”, “gaus”);
Double_t sb = 5.0; // sigma band for the fg integral (e.g. 3 or 5 or 7)

Int_t ch_min = 360, ch_max = 460; // “channel numbers”
TF1 *f1 = new TF1(“f1”, “gaus(0) + pol1(3)”, ch_min, ch_max);
f1->SetParameters(30000., 430., 10., 5000., 0.);
TFitResultPtr r1 = g1->Fit(f1, “WEMRS”);
TMatrixD c1 = r1->GetCovarianceMatrix();
TMatrixD c1g = c1.GetSub(0, 2, 0, 2);
Double_t s1, s1e;
fg->SetParameters(f1->GetParameters());
// fg->SetParErrors(f1->GetParErrors());
s1 = fg->Integral(fg->GetParameter(1) - sb * fg->GetParameter(2),
fg->GetParameter(1) + sb * fg->GetParameter(2));
s1e = fg->IntegralError(fg->GetParameter(1) - sb * fg->GetParameter(2),
fg->GetParameter(1) + sb * fg->GetParameter(2),
fg->GetParameters(),
c1g.GetMatrixArray());
std::cout << "gaus peak integral (mean +/- " << sb << " sigma) = "
<< s1 << " +/- " << s1e << std::endl;


You keep saying play the same game. I don’t get it. Of course, I know how to adapt this part to background as well.

The point is here when you say f1->SetParameters(30000., 430., 10., 5000., 0.);
As you know, as soon as we fix any of these parameters, we loose the error. I am asking what is your strategy of getting it back when they are still FIXED?

For instance, this part will not give me the error. I am avoiding to write all of the code. Because, I fixed the parameters where I store in my array of p_bg here.

double bcg_g1_error = bg->IntegralError( g1_X_LBE, g1_X_UBE , p_bg, cov_bg.GetMatrixArray()) / BinWidth;

The TF1::IntegralError method gets the parameter uncertainties (and their correlations) through its “covmat” input parameter.

I know this, but I think you do not get my question. I try again When you fix the parameters of any function, any function at all, is it possible to get the error in this method?

I know the answer: NO. Please, do correct me if I am wrong.

Here is a link you answered somebody else before: Integral Error in Signal - #6 by Wile_E_Coyote

Here is the adapted version of YOUR code below to get background:
As you know as soon as you say somewhere “fiz the parameter”, yo -u get ZERO error as normal.

{
60Codata.txt (17.3 KB)

const Bool_t use_partial_covariance = kFALSE; // kFALSE or kTRUE

const char *filename = “60Codata.txt”;

TGraph *g0 = new TGraph(filename, “%lg %lg %*lg”); // X axis = channels

g0->SetTitle(“g0;Channel number;Energy”); // energy calibration

TGraph *g1 = new TGraph(filename, “%lg %*lg %lg”); // X axis = channels

g1->SetTitle(“g1;Channel number;Counts”);

TGraph *g2 = new TGraph(filename, “%*lg %lg %lg”); // X axis = energies

g2->SetTitle(“g2;Energy;Counts”);

if ((g1->GetN() < 2) || (g2->GetN() < 2)) return; // just a precaution

TF1 *fg = new TF1(“fg”, “gaus”);

Double_t sb = 5.0; // sigma band for the fg integral (e.g. 3 or 5 or 7)

TF1 *bcg = new TF1(“bcg”, “pol1”);

Int_t ch_min = 360, ch_max = 460; // “channel numbers”

TF1 *f1 = new TF1(“f1”, “gaus(0) + pol1(3)”, ch_min, ch_max);

f1->SetParameters(30000., 430., 10., 5000., 0.);

f1->FixParameter(3, 5000); //for pol **************************************
f1->FixParameter(4, 0.); //for pol *******************************************

TFitResultPtr r1 = g1->Fit(f1, “WEMRS”);

TMatrixD c1 = r1->GetCovarianceMatrix();

TMatrixD c1g = c1.GetSub(0, 2, 0, 2);

TMatrixD c_bcg = c1.GetSub(3, 4, 3, 4);

Double_t s1, s1e;

fg->SetParameters(f1->GetParameters());

bcg->SetParameters(f1->GetParameters()+3);

// fg->SetParErrors(f1->GetParErrors());

s1 = bcg->Integral(fg->GetParameter(1) - sb * fg->GetParameter(2),

                fg->GetParameter(1) + sb * fg->GetParameter(2));

s1e = bcg->IntegralError(fg->GetParameter(1) - sb * fg->GetParameter(2),

                      fg->GetParameter(1) + sb * fg->GetParameter(2),

                      bcg->GetParameters(),

                      c_bcg.GetMatrixArray());

std::cout << "gaus peak integral (mean +/- " << sb << " sigma) = "

        << s1 << " +/- " << s1e << std::endl;

}

Additionally, try (just for “debugging”): c1.Print(); c1g.Print(); c_bcg.Print();

Have you run the code to see ZERO for the error? Not sure, how it will help. Of course, the matrix is full of ZERO.

You got also this post related to similar stuff as always. Apparently, it is a common issue for everbody.

None of these posts, as far as I noticed, are a solution to get integral error for the background when the parameters has to be fixed. Thus. another way to approach the problem is needed as I came to conclusion.

I got another idea of creating a histogram for the background to represent, then pull the error from the background histogram with IntegralAndError() method at the moment. The picture below is the correct representation of the background. if you give any flexiblity to the parameters, we make a distortion from this base line which makes it pol2 or something else indeed.

I saw this post from @moneta. I also try to understand that.

When fitting (make “background” parameters “free”), you can try to play with fit options: "EMRS", "PEMRS", "WEMRS", "LEMRS" ("L" only for histograms)

Ok. Now, we think the reason of getting different background estimations might be due to the fitting procedure we picked.

The TH1::Fit Method

Here are the results for background:

  1. Correct Result from histogram calculation:

Gross_Counts: 61171 ± 247.328 BCG_Counts: 15246 ± 693 Net_Counts: 45925 ± 735.812

  1. Result from fit option PEMRS:
    Integral range on X axis: xl: 400 xr: 462 g1_X_LowerBinEdge: 400 g1_X_UpperBinEdge: 463 Gross_area: 60885.7 ± 247.077 ( slightly UNDERESTIMATED)
    Bcg_area: 12807.1 ± 301.172 (UNDERESTIMATED) Net_area: 48078.6 ± 351.647 (OVERESTIMATED)

  2. Result from fit option EMRS:
    Integral range on X axis: xl: 400 xr: 462 g1_X_LowerBinEdge: 400 g1_X_UpperBinEdge: 463 Gross_area: 60885.7 ± 247.077 (SAME AS 2.)
    Bcg_area: 12807.1 ± 301.172 (SAME AS 2.)
    Net_area: 48078.6 ± 351.647 (SAME AS 2.)

  3. Result from fit option WEMRS:
    Bcg_area: 12422.8 ± 17.9087 (UNDERESTİMATED, it took longer time)

  4. Result from fit option LEMRS:
    Bcg_area: 12874.5 ± 296.152 (UNDERESTİMATED)

  5. Result from fit with fixed pol parameters and EMRS fit option:
    Bcg_area: 15062 ± 0 (similar as 1, the correct estimation)

  6. Geometric result: 15246 (from trapezium area of red area in my plot above.)

  7. My method: I calculate the coordiates of each end of the function limits and the slope of the background line in the code. Then, I create a function where I fix its parameters to those calculated results. Lastly, I use Integral() method without fitting anything. Namely, I read the integral over the function directly.

Results: Background integral: 15246 ± ??
(without error, EXACTLY SAME as 1st method.)

However, I don’t get the error because I haven’t fit anything. To make the IntegralError() method work, we need to fit and retrieve the covmatrix … from it. That’s why I got no ERROR retriaval.

NOTE: I have used same calculated coordinates for consistency to initialize the pol parameters throughout the methods 2-6 for fitting.

What is your comment on this?

Cheers.

In a comparison table, you should always add the area (and its error) of your fitted “signal”, too.

It is surprising to me that "EMRS" and "PEMRS" give the same results. I don’t know how you fill the histogram but, maybe you should apply (after filling, but before fitting): histo->ResetStats();

If you fill the histogram with weights different than 1, you should also try: "WLEMRS"

See the TH1::Fit method description for the available fit options.

I can do those. Bin width is 1. I have many functions to do everythiing separately.
I added Reset option, but it didn’t change the result for background. I can add gross area and net area, or counts as well.

I UPDATED first 2 results above in the previous one for you.

If you ask how I fill the gross area histogram, here is the function:


TH1D *CreateHisto(){

TH1D *hist = new TH1D(“histo”,“Gaussian Peak and Background”,1024,0,1024);
fstream file;

file.open("…Co60.txt", ios::in);

double Counts;
int channel=0;

while(!file.eof()){

file >> Counts ;
hist->Fill(channel,Counts);
//cout << channel << " " << Counts << endl;
channel++;
if(file.eof()) break;

}


Also, here is the part of an another function where I fill and create the background histogram.
I am not sure if SetBinContent and Fill methods make any difference.


TF1 *new_Bcg = new TF1(Form(“bg”),background ,Bcg_x1_LowE , Bcg_x2_UpE , 2); //pol1
new_Bcg->FixParameter(0, Bcg_a); //for pol
new_Bcg->FixParameter(1, Bcg_b); //for pol
cout << "\nBackground integral: " << new_Bcg->Integral( Bcg_x1_LowE , Bcg_x2_UpE) / BinWidth <<endl;

TH1D* bcg_histo = new TH1D(“bcg_histo”,“Background”,xr-xl,xl,xr);
//bcg_histo->Sumw2();

for (double i = xl; i < xr+1; i++){
bcg_histo->Fill(i,new_Bcg->Integral(i,i+1));
//bcg_histo->SetBinContent(i,new_Bcg->Integral(i,i+1));
//cout <<i << " " << new_Bcg->Integral(i,i+1) << endl;

}

bcg_histo->GetXaxis()->SetRange(xl,xr); // for fitting to focus on the limits

double BCG_Error ;
double BCG_counts = bcg_histo->IntegralAndError(xl, xr, BCG_Error , “”);
cout << "BCG_Counts: " << BCG_counts << " ± " << BCG_Error << endl;


After the “while” loop, add: hist->Sumw2(kFALSE); hist->ResetStats();

In principle, you should do the same for your bcg_histo (after the “for” loop): bcg_histo->Sumw2(kFALSE); bcg_histo->ResetStats();

As you fill the “bcg_histo” with “integral” values, for comparisons, you should probably add the "I" fit option to every fit.

I guess you mean after the end of the loops both in while and for. I can do that right now. You concern that it might have kept the previous result in the memory.

Now, it looks like these.

//////////////////////////////////////////////////////////
while(!file.eof()){
file >> Counts ;
hist->Fill(channel,Counts);
//cout << channel << " " << Counts << endl;
channel++;
if(file.eof()) break;

}
hist->Sumw2(kFALSE);
hist->ResetStats();
/////////////////////////////////////////////////////

for (double i = xl; i < xr+1; i++){
bcg_histo->Fill(i,new_Bcg->Integral(i,i+1));
//bcg_histo->SetBinContent(i,new_Bcg->Integral(i,i+1));
//cout <<i << " " << new_Bcg->Integral(i,i+1) << endl;

}

bcg_histo->GetXaxis()->SetRange(xl,xr);
bcg_histo->Sumw2(kFALSE);
bcg_histo->ResetStats();
////////////////////////////////////////////////////////////////////

I have only one fit. I don’t use any fit to create the bcg histogram. They are 2 different functions.

RESULTS ARE STILL THE SAME AS ABOVE. JUST 0.1 count different. It took much longer TIME with I option in the fit as I stated below.

BTW. Histogram integrals always use whole bins. So, for comparisions with function integrals, you should use:

Int_t xl = ..., xr = ...; // histogram bin numbers
... histogram->Integral(xl, xr) ...
Double_t xmin = histogram->GetBinLowEdge(xl); // left edge of the "xl" bin
Double_t xmax = histogram->GetBinLowEdge(xr + 1); // right edge of the "xr" bin
... function->Integral(xmin, xmax) ...

Don’t worry , I do it already. I am well aware of it. I attach the related part of the code here.
Could you please have a look?

Appreciated.
test1.C (11.0 KB)
Co60.txt (3.5 KB)