Gaussian fit + background in gamma spectroscopy

Newbie
1

Hello,
I am trying to perform a gaussian fit of a photopeak for several elements. I also want to include an exponential background noise but I have trouble in defining the final fit function. The Gaussian fit actually works fine but when I combine the two something weird shows up (see attached file).

ciao
ciao1918×916 107 KB

#include "TCanvas.h"
#include "TFile.h"
#include "TH1.h"
#include "TF1.h"
#include "TAxis.h"
#include "TStyle.h"
#include <string>
#include <iostream>

void analisidati_gamma1(
    std::string filename="Co60hq0.root"
    ,std::string histoname="hq0"
    ,std::string bkgfunc="expo"
    ,int xminGFit1=1120
    ,int xmaxGFit1=1240
    //,int xminGFit2=1250
    //,int xmaxGFit2=1400
    ){
    gStyle->SetOptFit(1111);

    TCanvas* can = new TCanvas("can","my canvas",800,800);

    // retrieve the histogram from the input file
    TFile* file = TFile::Open(filename.c_str());
    if(!file){
        std::cout << "File " << filename << " not found... check the input name!" << std::endl;
        return; 
    }
    TH1D* histo = (TH1D*) file->Get(histoname.c_str());
    if(!histo){
        std::cout << "histo " << histoname << " not found... check the input name!" << std::endl;
        return; 
    }
    
    
    // Get the TTree from the file
    TTree *nt = (TTree*) file->Get("nt");

    // Create a new histogram with 500 bins in the range [0, 10000]
    TH1D *histonew = new TH1D("histonew", "hq0new", 500, 0, 10000);

    // Draw the variable "qlong" from the TTree "nt" into the histogram "histonew"
    nt->Draw("qlong>>histonew", "ch==0");
    TH1D* histo_copy = (TH1D*)histonew->Clone();
    histo_copy->SetName("histo_copy");

    // Calibration
    TH1D *histo_cal = new TH1D("histo_cal", "Calibrated Histogram", 500, 0, 3000);
    nt->Draw("-32.4+qlong*0.1942>>histo_cal", "ch==0");

    // Gaussian fit

    const double normGaus=100; // Gauss normalization
    const double meanGaus=1700;    // Gauss mean
    const double sigmaGaus=500;

    // fit function
    std::string fitfunc = "gaus(0)";
    TF1* gaussFunc = new TF1("gaussFunc",fitfunc.c_str(),0,1000);

    gaussFunc->SetParName(0,"Gauss Norm.");
    gaussFunc->SetParName(1,"Gauss mean");
    gaussFunc->SetParName(2,"Gauss sigma");
    gaussFunc->SetParameter(100,normGaus);
    gaussFunc->SetParameter(1170,meanGaus);
    gaussFunc->SetParameter(500,sigmaGaus);

    //histo_cal->Fit(myfit,"","",xminGFit1,xmaxGFit1);

    //TFile* file_output = new TFile("file_output.root","recreate");
    //can->Write();
    //histo_cal->Write();


    // Define background exponential function
    TF1 *bkgFunc = new TF1("bkgFunc", "expo", 0, 3000);
    // Combine Gaussian and background functions
    TF1 *fitFunc = new TF1("fitFunc", "gaussFunc(2) + bkgFunc", xminGFit1, xmaxGFit1);
    // Set initial parameters for background exponential function
    bkgFunc->SetParameters(0, 1.0); // Exponential amplitude
    bkgFunc->SetParameters(1, -0.001); // Exponential decay constant
    // Set initial parameters for combined function
    fitFunc->SetParameters(normGaus, 1.0, -0.001);
     // Fit combined function to histogram
    histo_cal->Fit(fitFunc, "", "", 0, 3000);
return;
}
2

Hi Liz,

In order to test your script, you should also make the data file available.

By the way I do not see an exponential background. The stuff below channel
500 is just some detector pile-up and restrict your fit to channels > 500.

-Eddy

3

Hi, I’m sorry but unfortunately the file appears to be too large (I tried to zip it but it’s still to big).
Our professor told us to add an exponential background so that’s why we’re doing it.
Thank you for your reply

5

Hello,

There are a few problems with your code which can be easily modified.

  1. You try to fit a gaus + exp over the entire [0,3000] range this is why you have that red fitted line running across. When you define your fit function you limit it’s range:

but later, when you actually fit the histogram you overwrite this with [0,3000]:

If you check class reference for Fit you can find several fitting options including R which will limit the fitting range specified in the function range so down to your [xminGFit1, xmaxGFit1] range which you gave when you defined the fitFunc.

histo_cal->Fit(fitFunc, "R");

  1. I can see from these limits that you only try to fit the first peak with your gaus + exp which is not great when you have peaks this close together. You should at least use the [1000,1500] range and use two gaus functions for the 2 peaks + exp for the background.

  2. Your initial parameters for the Gauss peak are not good for two reasons.

  • You correctly use gaussFunc->SetParName(0, "Gauss Norm."); etc. to set the parameter names but when you try to set the values you don’t longer follow proper syntax.

This means that you try to set parameter number 100 to 100 and so on even though the gaus have 3 parameters, so you should use 0,1,2 like you did in the SetParName() lines. You didn’t notice but ROOT just ignored this and left your parameters on default 0,0,0 values and successfully found the actual values during fitting.

This line overwrites parameter number 0,1,2 with normGaus(=100), 1.0 and -0.001 and zeroes out the rest even though your 0,1 parameter belongs to exp and parameter number 2,3,4 belong to gaus. So you either go one by one using SetParameter() or use a single SetParameters() where you set every parameter in the correct order. Watch out, there is a difference between the two!

  • If you give better parameters to start with you have a better chance to get a proper fit. Your normGaus parameter is the amplitude and it is the height of the peak which you fit, so just estimate it reading the y-axis (e.g ~740 for first peak). Your meanGaus is the mean which you estimate by reading the x-axis position of your peak / you know it from a table. If you trust your calibration and know that these are Co60 peaks you can find in literature that peaks are at: 1173.2, 1332.5 keV energies (not 1700 which you gave). sigmaGaus parameter should be much lower, it corresponds to the width of your gaus. Play around with these and see the difference!

If you fix the range problem I mentioned in my first point you should get a normal fit.

Cheers,
LeX

6

Thank you so much for the advice, it turned out to be very useful!
Have a nice day :blush: