Calculating number of entries in a several peaks with different energies in a Histogram and calculating uncertainties and drawing error bars graph

[Mo922022]

Hello rooters,

I have 2 Root files from two different Simulation and as seen in the image is a comparison between energy spectrum of prompt gamma radiation produced from proton interactions. Since this is the last part of my thesis I want a help in this part please.

I want to Calculate the number of entries in some peaks in the histograms and since those entries are absolute number I want also to normalize the Histogram, moreover drawing error bar graph of uncertainties of both histograms and compare between them.
-peaks of interest are (0.71, 2.3, 4.43, 6.128) MeV.

I did a simple code to calculate number of entries in a specific peak but since each peak has an energy range I estimated binX1 and binX2 on X-axis by eye. but I want to do it by a scientific way to be accurate as possible.

• Here is the Code I have done
  DrawIntegral.cpp (2.8 KB)
• Here also the root files of Simulation
  Out_run_0.root (10.2 KB) Out_run_1.root (12.0 KB)

Any help would be appreciated
Thanks in Advance
Mo

[Axel]

You can use peak-finding algorithms that make this a bit more scientific. Or maybe better, fit the background, apply some cut on the signal height, and look for bins that have entries > background + minimal_signal_height. Then fit these ranges (Gaussians?); take the fit integral and that’s the number of entries in the peaks.

You can normalize a histogram by calling hist->Scale(1./hist->GetSumOfWeights()) (which excludes entries in under- and overflow bins).

[Mo922022]

Thanks a lot for your reply

I tried the peak Finding algorithms with TSpectrum class if u mean this method and this does not help too much or I could say , it is only finding the peak and I already know the peaks I want. I tried also fitting by TF1 and it somehow works but still I have to estimate range of Xmin and Xmax by eye here is the point to be more scientific is there is a possibility to calculate peak integral by calculating FWTM or some thing like this and this could be helpful.

In the following code u will find the code I have done to generate the following Canvas.

DrawPeakIntegral.cpp (7.0 KB)

in the Following lines you will see calculation of FWHM or FWTM but the problem it gives for whole Spectrum not for a specific peak.

double xmin = hEnergyGamma->GetMinimum(4.3);
double xmax = hEnergyGamma->GetMaximum(4.6);
int bin1 = hEnergyGamma->FindFirstBinAbove(hEnergyGamma->GetMaximum()/2); // divide by 10 for FWTM
int bin2 = hEnergyGamma->FindLastBinAbove(hEnergyGamma->GetMaximum()/2); // divide by 10 for FWTM
double FWHM = hEnergyGamma->GetBinCenter(bin2)- hEnergyGamma->GetBinCenter(bin1);
std::cout << "#####$$$$$ FWHM = " << FWHM << std::endl;

[Axel]

Once you fit you can use the fit’s TF1::Integral(). For fitting you indeed need to give starting values for the parameters; for any “reasonable” value you should arrive at more or less the same fit result. Both elements are bringing the “scientific stability” you’re after

[Mo922022]

Okay, Thanks a lot