fitting a sloped gaussian

35 vues (au cours des 30 derniers jours)
Scarlet Passer
Scarlet Passer le 17 Fév 2021
Commenté : dpb le 18 Fév 2021
Hello!
So I have some data from a compton scattering experiment and I want to fit my peak with a gaussian and find the standard deviation, the amplitute, and the area under the curve. I've tried the 'gauss3' fit and it works pretty well but I'm not sure where to go from there. Because the background is sloped I don't know how to tackle the area without getting unwanted space from beneath the sloped background. All I want is the area of the curve itself. This is what it looks like:
Thanks in advance!
  5 commentaires
Alex Sha
Alex Sha le 18 Fév 2021
If use the combination of 5 Gauss function, the result will be much better:
y = a1*exp(-((x-b1)/c1)^2)+a2*exp(-((x-b2)/c2)^2)+a3*exp(-((x-b3)/c3)^2)+a4*exp(-((x-b4)/c4)^2)+a5*exp(-((x-b5)/c5)^2)
Root of Mean Square Error (RMSE): 10.2702586328186
Sum of Squared Residual: 108115.167694609
Correlation Coef. (R): 0.995045924385973
R-Square: 0.990116391637135
Parameter Best Estimate
---------- -------------
a1 92.7498096695778
a2 -1769.94724857124
a3 -145.468199223536
a4 60.6523446676163
a5 151714694.184991
b1 199.060149101054
b2 -61.3707462724171
b3 46.4480539913997
b4 456.40444211691
b5 -2441.31620849057
c1 17.1540546396503
c2 95.5557434478509
c3 26.7272034725448
c4 -31.7643462462293
c5 -711.910806560653
dpb
dpb le 18 Fév 2021
But there are "issues" here...
a2, a3 < 0
b2, b5 < 0
c4, c5 < 0
and no background term estimated at all, either. I didn't have time to do anything more last night, sorry...

Connectez-vous pour commenter.

Réponses (1)

dpb
dpb le 18 Fév 2021
Modifié(e) : dpb le 18 Fév 2021
First pass...overall not so great; need to do more to figure out how to model background. My spectroscopy background has dealt only with gamma-spec so not familiar with what the accepted methods are for Compton scattering data; the decaying exponential can produce a nice fit overall, but the baseline estimate is going to be very variable from run to run I suspect. I tried a couple ranges; your limited LH on the peak and another wider with similar results. But, for starters,
N1=find(energy>400,1); % this returns the 464 you used, btw...by happenstance
E=energy(N1:end); S=int50(N1:end); % shorthand variables for the chosen range
% define fitting function including exponential decay background
fnPk=@(a,b,a1,b1,c1,x) a*exp(-b*(x-E(1)))+a1*exp(-((x-b1)/c1).^2);
fnBG=@(a,b,x) a*exp(-b*(x-E(1))); % shorthand for the background term for later...
f50=fit(E,S,fnPk,'Lower',[25 0 0 0 0],'StartPoint',[50 0.001 60 450 20]);
This returns
>> f50
f50 =
General model:
f50(x) = a*exp(-b*(x-E(1)))+a1*exp(-((x-b1)/c1).^2)
Coefficients (with 95% confidence bounds):
a = 25 (23.31, 26.69)
b = 0.004602 (0.004254, 0.00495)
a1 = 54.7 (52.99, 56.41)
b1 = 456.1 (455.5, 456.7)
c1 = 24.55 (23.51, 25.59)
>> figure
>> plot(E,S), hold on, plot(f50)
>> hBG=plot(E,fnBG(fitPK.a,fitPK.b,E),'k-');
>> hLG=legend; hLG.String={'Fitted model','Background'};
results in
A wider swath using N1=find(E>350,1) resulted in very similar fitting; a little higher background level but still significantly below the visual location when estimated in the combined model.
Had a similar issue years and years ago in another world with a NaI scintillation setup where were trying to use a quadratic for background; we ended up fitting the two independently with OLS for the background quadratic and the gaussian separately with a combined Marquardt routine. It improved things some, but the background estimation was always problematical for reproducibility from repitition to repition of measurements of the same sample.
That was for an online industrial instrument; we ended up finally having to average successive measurments as simply didn't have sufficient computing power at the time to do anything more sophisticated. This, of course, was 30 years ago by now...
Good luck, I'll do a very quick search for Compton spectra and see if anything shows up if get a chance...
  1 commentaire
dpb
dpb le 18 Fév 2021
The wider before I delete it..

Connectez-vous pour commenter.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by