Gauss and Pol1 Fitting Results

Hi guys,
Recently I started to doubt myself, so I need help.

Imagine you have gamma spectrum with one isolated single peak like 137Cs peak at 661keV.
Anyways. If we want to get the integral and its error of the peak, how to do it properly? In the end, I will compare the results to getting the answer from histogram bin contents directly using either Covell method or Total Peak Area method.

The fit section of the code:

I tried to simplify my question by making up the limits like a,b,c,d,e,f. I also skip other lines which is not the focus of my question.

TF1 *g1 = new TF1(“g1”,“gaus”, a, b );
TF1 *pol1 = new TF1(“pol1”,“pol1”, c, d );
TF1 *total = new TF1(“total”,“gaus(0)+pol1(3)”, e , f );
Then we fit this onto a histogram of the data.

h1->Fit(g1,“R”);
h1->Fit(pol1,“R+”);
h1->Fit(total,“R+”);

cout<<“g1 from integral : " << g1->Integral( p1-3p2, p1+3p2 ) <<endl;
cout<<” pol1 from integral: " << pol1->Integral( p1-3p2, p1+3p2 ) <<endl;
cout<<" total from integral : " << total->Integral( p1-3p2, p1+3p2 ) <<endl;

How is that possible that I got the integral of total bigger number than the Gauss integral here?

Also, why do we define the total function like “Gauss+pol” instead of “Gauss-pol” because we need to deduce the background?

Basically, what am I missing?

Cheers.

You didn’t put any numbers in your example, so I don’t know what you are seeing, but it is usually obvious that a sum model sees more than it’s components.
Imagine you have two processes, signal and background. The background process generates 100 events, the signal process generates 10.
In total, you will observe 110 events.

If you fit this, you need a signal + background model, because you have both signal and background in the data. Now somebody asks how much signal you observe. So you fit the combined model to the (combined) data, and then you extract how much the signal model sees.
Results will be:

  • Combined: 110
  • Background: 100
  • Signal: 10

Obviously with some fluctuations within uncertainties.



Hi,
I put the number in txt mode as attached.
CS137.txt (10.0 KB)
It’s a point source data collected by NaI detector, that’ it. Simple one. No need to go through the possible background sources and so on. My point is only to assume the background as flat polynomial in 1st order and deduce it from the range of interest of the peak. Then find the net count. By doing these, I use bins and their contents due to the method I choose in uncalibrated data in terms of energy. Since I know which peak is which one, that’s fine and more useful.

The part for fitting one Gaussian.

void FittingOneGussianFunctions(TH1F *h1, int LL1, int UL1, double & p2){

////////////////////Defining Functions //////////////////////////////////////////////////////////////
//The limits are put t get exact function shapes, so they might be narrow.
int flexibility1=20;
TF1 *g1 = new TF1(“g1”,“gaus”, LL1, UL1 );
TF1 pol1 = new TF1(“pol1”,“pol1”, LL1-flexibility1, UL1+flexibility1 ); //POLYNOMIAL BACKGROUND
TF1 total = new TF1(“total”,“gaus(0)+pol1(3)”, LL1 , UL1 );//TOTAL FUNCTION
/////////////////////////////////////////////////////////////////////////////////////////////////////////
//TH1F h = (TH1F)h1->Clone(“h”);
h1->SetLineColor(kBlack);
g1->SetLineColor(kBlue);
pol1->SetLineColor(kBlack);
total->SetLineColor(kRed);
h1->Fit(g1,“R”);
h1->Fit(pol1,“R+”);
h1->SetStats(000000000);
///////////////////////////////////////////////////////////
// To store the parameter values in to an ARRAY called “par[5]”
double par[5];
g1->GetParameters(&par[0]);
pol1->GetParameters(&par[3]);
total->SetParameters(par); //first get the parameter from g1 and pol1 , then set them for total function as parameters.
h1->Fit(total,“R+”); //Now, I am able to fit and draw this function since I assigned its values which were taken from the first gaussian and pol1 fits
//double chi2 = fitting->GetChisquare();
// value of the first parameter
double p0 = g1->GetParameter(0); //height of gaussian
double p1 = g1->GetParameter(1); //mean value of gaussian
p2 = g1->GetParameter(2); //if this is sigma, width of gaussian=8
sigma
cout<<"height value from fit= "<<p0 << endl;
cout<<"mean value from fit= " <<p1 <<endl;
cout<<"sigma value from fit= " <<p2 <<endl;
cout<<"FWHM from fit= "<<2.35
p2 <<endl;
////////////////
double p3 = pol1->GetParameter(0);
double p4 = pol1->GetParameter(1);

// error of the first parameter
double e0 = g1->GetParError(0);
double e1 = g1->GetParError(1);
double e2 = g1->GetParError(2);
//////////////////////////////
double e3 = pol1->GetParError(0);
double e4 = pol1->GetParError(1);
///////////////////////////////////////////////
//For interested gaussian peak ONLY
double* gaussParameters = new double[3] ;
gaussParameters[0] = p0;
gaussParameters[1] = p1;
gaussParameters[2] = p2;
double* MyConvarienceMatrix =new double[9];
MyConvarienceMatrix[0]=e0 ;
MyConvarienceMatrix[1]=0 ;
MyConvarienceMatrix[2]=0 ;
MyConvarienceMatrix[3]=0 ;
MyConvarienceMatrix[4]=e1 ;
MyConvarienceMatrix[5]=0 ;
MyConvarienceMatrix[6]=0 ;
MyConvarienceMatrix[7]=0 ;
MyConvarienceMatrix[8]=e2 ;
/////////////////////////////

int HistoBinning=1;
double MyIntegral = g1->Integral( p1-3p2, p1+3p2 ) / HistoBinning ; //p4 is mean value and p5 is sigma value. Gauss range is about 8sigma.
double IntegralError = g1-> IntegralError(p1-3
p2 , p1 + 3p2 , gaussParameters, MyConvarienceMatrix) / HistoBinning ;
cout<<"Fitting g1 integral and error: " << g1->Integral( p1-3
p2, p1+3p2 ) <<"\t"<<IntegralError <<endl;
cout<<"Fitting pol1 integral and error: " << pol1->Integral( p1-3
p2, p1+3p2 ) <<endl;
cout<<"Fitting total integral and error: " << total->Integral( p1-3
p2, p1+3*p2 ) <<endl;
//SET gaussian ranges for drawing
g1->SetRange(LL1-flexibility1, UL1+flexibility1);
pol1->SetRange(LL1-flexibility1, UL1+flexibility1);
total->SetRange(LL1-flexibility1, UL1+flexibility1);
////////////////////////////////////////
g1->Draw(“same”);
pol1->Draw(“same”);
total->Draw(“same”);

TLegend *legend1 ;
Double_t xm1=0.6, ym1=0.6, xm2=0.8 , ym2=0.88;
legend1 = new TLegend(xm1,ym1,xm2,ym2);
legend1->SetHeader("^{137}Cs 661.7 keV peak");
legend1->AddEntry(g1,“Gross Area”,“l”);
legend1->AddEntry(pol1,“Background Area”,“l”);
legend1->AddEntry(total,“Net Area”,“l”);
legend1->SetBorderSize(0);
legend1->SetTextSize(0.03);
legend1->Draw();

}

In your 2nd code, what exactly this part :" dx = (g2->GetX()[ch_max] - g2->GetX()[ch_min]) / (ch_max - ch_min);" adjusting? It looks like the slope of energy vs channel plot. First integral you find is divided by the dx, then the numbers are too different. You go 5sigma from the mean on both side to get the integral, do you take the portion in to account to retrieve the whole count. Let’s say if you used 1sigma, it would have covered 68.3% of the actual data.

In the first code, integral is taken from gauss+background, but in the second one it is taken from gauss only. My question is still the same: why we dont take the integral from gauss-background?

Fits look perfect, by the way.

{
  TTree *t = new TTree("t", "t");
  t->ReadFile("CS137.txt", "counts/D");
  t->SetEstimate(-1);
  t->Draw("counts:Entry$", "", "L");
#if 1 /* 0 or 1 */
  TGraph *g = (TGraph*)gPad->GetPrimitive("Graph"); // 2D
  TF1 *f = new TF1("f", "gausn(0) + pol1(3)", 210., 280.);
  f->SetParameters((12000. * 8. * 2.5), 245., 8., 100., 0.);
  g->Fit(f, "WEMR");
  std::cout << "gaus peak integral = " << f->GetParameter(0)
            << " +/- " << f->GetParError(0)
            << std::endl;
#endif /* 0 or 1 */
}

I saw what you tried to do earlier. You want to compare gauss integral and adjusted gauss+background integral. Which one is correct or which one is preferable or why we need the same one if only the gauss integral works?
Maybe, I missed a point.

My codes are so primitive and naive. You guys are so professional. I guess,sadly, it is my lack of computer skill dated back to my undergraduate years.

In the code above, how did you guess quickly the set parameters, and how does the first multiplication part as an initialization parameter become the integral of the fitting after fit ?

Can anyone tell me which one is actually corresponding to the net count area?

TF1 *g1 = new TF1(“g1”,“gaus”, LL1, UL1 );
TF1 *pol1 = new TF1(“pol1”,“pol1”, LL1, UL1 ) ; //POLYNOMIAL BACKGROUND
TF1 *total = new TF1(“total”,“gaus(0)+pol1(3)”, LL1 , UL1 );//TOTAL FUNCTION

Fitting g1 integral : 207802
Fitting pol1 integral: 4452.39
Fitting total integral: 209342

It seems from the result of integrals for the same bin range, total function is really the the total of the two, which then it means nothing to me. Net count is the gross area,g1, minus background, pol1. Right?
It sounds like I draw total function first and then background and then gaussian. Then, I kinda try to rretrieve gaussian and background from total function one by one. ıN MY MİND, it should be like we draw gaus and background, then we get total function, which is the deduction of them.

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