Welcome!

This is a Mathematica notebook. Mathematica is a powerful symbolic and computational package, which you can download to your own computer from the CMU Licensd Software site. You’ll find text like this interwoven with blocks of code in a different font. You can execute a block of code by putting your cursor in the block and hitting shift-enter. The results - numerical, analytical, graphical, whatever - will show up below. Try it out by excuting this code:

Manipulate[ParametricPlot[{Sin[a x],Sin[b x]},{x,0,2π}],{a,1,5},{b,1,5}]

What we’re doing

The theory of hyperfine splitting says that Mn has electronic transition energies hν = ΔE = 2(Subscript[μ, B])B + A(Subscript[m, I]) where (Subscript[m, I]) is the nuclear spin quantum number. Therefore, to a first approximation, we would expect to see see six equally spaced resonance peaks. Even though the data aren’t perfect, we’ll also assume that each transition absorbs the same amount of the driving microwaves, leading to equal peak amplitudes, and we’ll approximate equal lifetimes, leading to a uniform linewidth for all the peaks.

Your ESR spectra show you the derivative of the resonance curve, so first we’ll create a function from the derivative of six Lorenzians using the assumptions above.

Next, we’ll import your data.

Next, we’ll perform a nonlinear fit of our Lorenzian function to your data, and plot the fit and residuals until we get it right.

Next, we’ll do some statistics to figure out how good our fit is and how certain we are of the fitting parameters (particularly the linewidth, which we’re after).

Setup and Import

This creates a six-peak Lorentzian derivative function P6 of x, with the first peak centered at (Subscript[x, o]) and subsequent peaks spaced uniformly by δ after the first peak, with uniform peak amplitude A and width Γ, and a constant offset b for the whole function.

P6[x_,x0_,δ_,A_,Γ_,b_]:=(A D[PDF[CauchyDistribution[x0,Γ], x],x ]+A D[PDF[CauchyDistribution[x0+δ,Γ],x],x]+A D[PDF[CauchyDistribution[x0+2δ,Γ], x],x]+A D[PDF[CauchyDistribution[x0+3δ,Γ], x],x]+A D[PDF[CauchyDistribution[x0+4δ,Γ], x],x]+A D[PDF[CauchyDistribution[x0+5δ,Γ], x],x])+b
P6[x,x0,δ,A,Γ,b]//Simplify

Export your oscilloscope’s collected data as a csv (see the sign on the wall). Make sure you have just two columns in your spreadsheet, then put the path to it here. Check that output looks like your data.

Data=Import["C:\\Users\\Owner\\Documents\\Spring 2011\\Modern Physics Laboratory\\ESR\\FILE.csv"]//TableForm

Fit the Data

This fits your data to the six-peak function. You’ll need to provide starting guesses (REPLACE_ME) for the fit parameters. The output gives you a list of the fit parameters, and Mathematica’s “error bars” on those parameters. We could take those at face value as our uncertainties... but what are they? Where do they come from?

DataFit=NonlinearModelFit[Data,P6[x,x0,δ,A,Γ,b],{{x0,REPLACE_ME},{δ,REPLACE_ME},{A,REPLACE_ME},{Γ,REPLACE_ME},{b,REPLACE_ME}},{x}]
DataFit[{"BestFitParameters","ParameterErrors"}]

Now let’s plot! The fit is in red and your data is in blue. If the fit follows your data well, great! If it’s totally wrong, then play with the starting values above and redo the fit until you get it. (P.S. You might need to adjust the plot ranges below, depending on how big your gain was.)

Show[ListPlot[Data],Plot[Evaluate[P6[x,x0,δ,A,Γ,b]/.DataFit["BestFitParameters"]],{x,-11,50},PlotStyle->Red,PlotRange->{-.6,.6}],PlotRange->{-.6,.6}]

Now we plot the residuals. In principle, these are supposed to be random....

ListPlot[Transpose[{Data[[All,1]],DataFit["FitResiduals"]}],PlotRange->Full]

Do some statistics

The residuals are not random, but we still claim the hyperfine splitting theory is a good approximation and will continue to use this model to represent our data.

We need to find uncertainties on our fit parameters somehow. We could try to wring these out of Mathematica, from the “error bars” you found above, but we don’t know what exactly those represent. We’ll have to find them by hand.

Unfortunately, most of the techniques we use assume that the errors in the data are random. If your residuals above looked random, good for you, but if they were anything like ours, you’ll have to work around it somehow.

Let’s start with the basics. (χ^2) is a quantity that we use to assess goodness of fit. We take the summed amount that the data differs from the fit (and square it, because variation in opposite directions should add, not cancel), and normalize it by the variation that we expect in each point due to random error (σ). A related quantity is reduced (χ^2), or (Subscript[(χ^2), r]), where we divide through by the number of “degrees of freedom” - essentially, data points minus fit parameters. Since bigger data will intrinsically have more total random error, and more fit parameters will tend to produce a tighter fit, (Subscript[(χ^2), r]) gives us a convenient “universal” measure of goodness of fit. We generally expect (Subscript[(χ^2), r]) to be around 1. Too much lower and something’s suspicious (often an artificially high σ). Too much higher and we’re tempted to say the model doesn’t fit the data.

So let’s compute (χ^2).  First we sum the squared residuals from our fit above. Then we divide by a baseline uncertainty (σ) for each data point (one way you can get this is by looking at the fluctuations in the tail of your spectrum, where you expect the signal to be constant). Then, to get (Subscript[(χ^2), r]) we divide by the degrees of freedom. (Don’t be distressed if your numbers are a little high - remember that your errors aren’t necessarily random.)

σ = REPLACE_ME;
χsq = Total[((DataFit["FitResiduals"])^2)]/(σ^2)
χsqreduced=χsq/(Length[Data]-Length[DataFit["BestFitParameters"]])

To find uncertainties on our fit parameters, we think of (χ^2) as a function of our fit parameters. (In this case, we have four parameters, so a function on a four-dimensional space.) A best fit finds the global minimum of (χ^2) (in principle). If the local curvature near the minimum is very high, it means that moving a parameter just a little will kick (χ^2) way up, so we can be pretty confident in the value of the parameter. On the other hand, if the curvature is small, we can move the parameter around a lot without changing the goodness of fit too much, which gives us low confidence in the parameter. We quantify this by saying that the uncertainty on our parameter is the amount we need to change it to move (χ^2) up by 1. (Remember, (χ^2), not (Subscript[(χ^2), r])!) Remember that this assumes random errors, so the errors on our parameters should be uncorrelated: in principle, no matter which direction we move in parameter-space, we should hit our increased (χ^2) at the same distance out on some sort of n-dimensional circle.

Our errors aren’t random, though. To try to capture the variation, we’re going to pull a trick that’s not entirely justified. We’ll move out from our best-fit values of Γ and δ, but then re-fit the other parameters and re-calculate (χ^2), and play this game until (χ^2) increases by 1. You can use the code below to calculate fits for particular fixed values of the parameters, and the code above to re-calculate (χ^2).

DataFitΓ=NonlinearModelFit[Data,P6[x,x0,δ,A,CHANGE_ME,b],{{x0,REPLACE_ME},{δ,REPLACE_ME},{A,REPLACE_ME},{b,REPLACE_ME}},{x}]
Total[((DataFitΓ["FitResiduals"])^2)]/(σ^2)

DataFitδ=NonlinearModelFit[Data,P6[x,x0,CHANGE_ME,A,Γ,b],{{x0,REPLACE_ME},{δ,REPLACE_ME},{A,REPLACE_ME},{b,REPLACE_ME}},{x}]
Total[((DataFitδ["FitResiduals"])^2)]/(σ^2)

Next Steps

Do you see any pattern in your residuals that you think can be fixed by adding another fitting parameter? (For instance, the assumption of uniform peak height or spacing could have caused non-random distortion.) If you can re-fit your data with extra parameters, can you think of physical justifications for those parameters?