#!/usr/bin/env python # ----------------------------------------------------------------------------------- # # # ROOT macro for reading text (ASCII) files and fitting data written. # # This exercise is about more advanced fitting of data. The outset is the simplest # possible function that even remotely does the job, which is then expanded to # accommodate all the features present in the data fitted. # The data is from a harmonic oscillator, in this case a spring hanging with # a weight, where the air drag gives a light damping. # # The formula for the position (d) as a function of time (t) of a simple (underdamped) # harmonic oscillator can be found here: d(t) = A * sin(w*t+phase) * exp(-gamma*t) # # Author: Troels C. Petersen (NBI) # Email: petersen@nbi.dk # Date: 9th of December 2015 # # ----------------------------------------------------------------------------------- # from ROOT import * from array import array # Setting what to be shown in statistics box: gStyle.SetOptStat("emr") gStyle.SetOptFit(1111) SavePlots = False verbose = True Nverbose = 10 # ----------------------------------------------------------------------------------- # # Define files, lists, counters and read data: # ----------------------------------------------------------------------------------- # file = "data_HarmOsc1.txt" time0 = 0.0 # For ensuring that time starts at 0.0s! time = [] dist = [] Nread = 0 # Loop over files and open them with open( file, 'r' ) as f : # Loop over lines and extract variables of interest for line in f: line = line.strip() columns = line.split() if (len(columns) == 2) : # Put the values read into lists (defining start at t=0): if (Nread == 0) : time0 = float(columns[0]) time.append(float(columns[0]) - time0) dist.append(float(columns[1])) if (verbose and Nread < Nverbose) : print " Time = %6.3f Dist = %6.3f "%(time[-1], dist[-1]) Nread += 1 else : print " Error: Strange values/lines read from data file!" f.close() print " Number of entries read: %d Time of last read: %6.3f"%(Nread, time[Nread-1]) # ----------------------------------------------------------------------------------- # # Plot the data: # ----------------------------------------------------------------------------------- # t = array( 'f', [0.0]*Nread ) et = array( 'f', [0.0]*Nread ) d = array( 'f', [0.0]*Nread ) ed = array( 'f', [0.0]*Nread ) # Fill the arrays from the lists and also assign uncertainties: # ------------------------------------------------------------- for i in range( Nread ) : t[i] = time[i] et[i] = 0.0 d[i] = dist[i] ed[i] = 0.042 # The size of the error on each point - initially WRONG! if not ((-0.001 < time[i] < 100.0) and (-5.0 < dist[i] < 5.0)) : print "Warning: Strange value for time and/or dist!", i, time[i], dist[i] # Make a graph of the data: graph = TGraphErrors(Nread, t, d, et, ed) # Make a white canvas and draw the data in: c0 = TCanvas("c0","", 20, 20, 1400, 700) graph.GetXaxis().SetRangeUser(0.0, time[Nread-1]+0.25) # Set the plotting range! graph.SetLineColor(kBlue) graph.SetLineWidth(2) graph.SetMarkerStyle(20) graph.SetMarkerColor(kBlue) graph.SetMarkerSize(0.5) graph.Draw("AP") c0.Update() if (SavePlots) : c0.SaveAs("Fit_HarmOsc1_v1.pdf") raw_input( ' ... Press enter to exit ... ' ) # ----------------------------------------------------------------------------------- # # # Looking at the plot of data from the file "data_HarmOsc1.txt", there is little doubt # that it is oscillating, and looking closer you will also see the slight damping. # # Questions: # ---------- # 1) Look at the data file and plot and see if you can by eye (or simple fits) estimate # the size of the uncertainty of the points. It is not easy, but you should be able # to get it to within a factor of 2-3. Try for 5-10 minutes and discuss it with your # neighbor before reading on! # # - - - - - - 5-10 minutes (success or failure) later - - - - - - # # If you didn't know how to estimate this uncertainty, then try to zoom in on a very # small part of the curve, where it should be possible to fit it with a line, or # possibly a 2nd or 3rd degree polynomial. Do so, and since you know that for a small # enough range of the data, this will be a reasonable PDF to use, extract the error # from the RMS of the residuals, which is roughly equivalent to choosing the errors, # that gives you a reasonable Chi2. If you want to check this, then using this error # for all points, see if it still fits well in some other (small) range. # # 2) Once you have tried, set the error to 0.0035, and try to fit a damped harmonic # oscillator. You will have to write the function yourself. Do so and run the fit # before reading on. # # - - - - - - 5-10 minutes (success or failure) later - - - - - - # # Did you manage to write the fitting function? Also, did you remember to put in a # parameter taking care of the offset from zero? If not, then look at the bottom of # this page, and "borrow" the parts that you need. Now run this fit, and see what # happens. # # - - - - - - 5-10 minutes (success or failure) later - - - - - - # # 3) Did the fit converge? I imagine that it didn't (thought it might have), and my first # guess on why (apart from copying the below wrongly) would be initial parameters. You # need to set the initial parameters quite accurately, for the fit to work. Think about # this, i.e. in a 5D parameter space, how hard it is to find just the correct frequency # when all you got is a Chi2 value, and being just 5% wrong or so will give you nothing # of value? # So now try to evaluate what some good initial values for the fit would be. You can # start by simply making educated guesses, but if that fails, you can instead draw the # function choosing some parameters, until it starts looking like the data you have. # Note that especially the frequency requires a good initial value. # # 4) Try to fit only the range 20s-22s, and see how significant the damping is. Would # you from only two seconds data be able to tell, that there is a damping, and if # so, by how many sigma? # # # # ---------------- Now consider the data file "data_HarmOsc2.txt" --------------- # # 5) Now consider the data file "data_HarmOsc2.txt", and set the uncertainty on the # distance to be 0.0025. Again, make a simple harmonic oscillator fit run. # Now your job is to expand on the fitting function and introduce terms to include # various effects and thus reduce the Chi2. # Set your fit to the range [0.005,36.005], and see how low a Chi2 you can get. # # 6) Can you detect a change in the distance uncertainty as a function of time? # I.e. does the errors needed to get a reasonable/constant Chi2 get better or worse # if you fit small ranges of time early and late in the experiment. # # # # ---------------- Now consider the data file "data_HarmOsc3.txt" --------------- # # 7) The data file "data_HarmOsc3.txt" contains a different type of damping in the # oscillation. Fit this, and determine at which point the oscillation stops. # # # ----------------------------------------------------------------------------------- # # # You need to write this functional form of a damped harmonic oscillator: # def func_osc(x, par) : # return par[0] + par[1]*sin(par[2]*x[0]+par[3])*exp(-par[4]*x[0]) # # fit = TF1("fit", func_osc, 0.0, time[Nread-1], 5) # fit.SetLineColor(kRed) # fit.SetLineWidth(2) # fit.SetNpx(2000) # This is just to make sure that the function is drawn smoothly enough! # graph.Fit("fit", "R")