#!/usr/bin/env python # ----------------------------------------------------------------------------------- // # # ROOT macro for introducing ROOT with examples. # # Author: Troels C. Petersen (NBI) # Date: 7th of October 2015 (latest update) # # Run by: ./RootIntro.py # # ---------------------------------------------------------------------------------- // # First, import the modules you want to use (here done differently than in PythonIntro.py): from ROOT import * from array import array # Arrays are much like python lists, but ROOT requires arrays not lists! import math # Next, we inform ROOT what plotting style it should use for statistics and fitting box: gStyle.SetOptStat("emr") # Options for statistics box: "e" = entries, "m" = mean, "r" = RMS gStyle.SetOptFit(1111) # Options for fitting box. r = TRandom3() # Random generator r.SetSeed(1) # Set a random seed (but a fixed one - more on that later.) SavePlots = False # For now, don't save plots (once you trust your code, switch on) verbose = True # For now, print a lot of output (once you trust your code, switch off) Nverbose = 10 # But only print a lot for the first 10 random numbers #---------------------------------------------------------------------------------- # Histogram with fit (1D): #---------------------------------------------------------------------------------- # Statistics and fitting results replaced in (can be omitted): gStyle.SetStatX(0.92); # Top right corner. gStyle.SetStatY(0.92); # Define histograms (name, title, number of bins, minimum, maximum): # It is a good idea to let the name be the same as defined outside. Hist_x = TH1F("Hist_x", "Hist_x", 100, -5.0, 5.0) Hist_y = TH1F("Hist_y", "Hist_y", 100, -5.0, 5.0) # Loop to get some random values and fill them into histogram: Npoints = 10000 # Number of random points produced for iexp in range( Npoints ) : x = r.Gaus()*0.8-0.5 Hist_x.Fill(x) y = r.Gaus()*1.3+0.5 Hist_y.Fill(y) if (verbose and iexp < Nverbose) : print " Histogram input value %2d: x = %6.3f y = %6.3f"%(iexp, x, y) # Plot result: # ------------ canvas = TCanvas( "canvas", "canvas", 50, 50, 1200, 600 ) Hist_x.SetTitle("Distribution of Gaussian numbers") Hist_x.GetXaxis().SetRangeUser(-5.0, 5.0) Hist_x.GetXaxis().SetTitle("Random numbers") Hist_x.GetYaxis().SetTitle("Frequency") Hist_x.SetLineColor(kBlue) Hist_x.SetLineWidth(2) Hist_x.Draw("e") # The option "e" shows errors (Poisson!) # Fitting histogram (with predefined function): fit_x = TF1("fit_x", "gaus", -2.0, 1.0) fit_x.SetLineColor(kBlue-8) fit_x.SetLineWidth(3) Hist_x.Fit("fit_x", "") # Drawing a second histogram: Hist_y.SetLineColor(kRed) Hist_y.SetLineWidth(2) Hist_y.Draw("same") # The option "same" makes it plot on top # Example of how to get e.g. means from a histogram and result of fit: print "Means: mu_x = %6.3f mu_y = %6.3f"%(Hist_x.GetMean(), Hist_y.GetMean()) print "Fitted mean of x: mu_hat = %6.3f +- %5.3f"%(fit_x.GetParameter(1), fit_x.GetParError(1)) # Legend: leg = TLegend( 0.15, 0.70, 0.38, 0.85 ) leg.SetFillColor(kWhite) leg.SetLineColor(kWhite) leg.AddEntry(Hist_x, " Gaussian (#mu = -0.5)", "L") leg.AddEntry(fit_x, " Fit with Gaussian to x", "L") leg.AddEntry(Hist_y, " Gaussian (#mu = +0.5)", "L") leg.Draw() canvas.Update() if (SavePlots): canvas.SaveAs("Histogram.pdf") #---------------------------------------------------------------------------------- # Graph with fit (1D): #---------------------------------------------------------------------------------- # Statistics and fitting results replaced in: gStyle.SetStatX(0.52); # Top left corner. gStyle.SetStatY(0.86); # Define points for a graph with errors: # Graphs in ROOT requires arrays, which is why these were included at the top! Ndata = 51 x = array( 'f', [0.0]*Ndata ) y = array( 'f', [0.0]*Ndata ) ex = array( 'f', [0.0]*Ndata ) ey = array( 'f', [0.0]*Ndata ) # Having defined four arrays filled with zeros, below we fill it with values: for i in range ( Ndata ) : x[i] = -1.25 + 0.05*i y[i] = 2.0*x[i]*x[i]*x[i] + x[i]*sin(x[i]) - 2*x[i] - 2.0 if (x[i] > 0.0) : y[i] += 0.5 # Put in a discontinuity by hand! y[i] += r.Gaus(0.0, 0.1) # Adding a Gaussian variation of 0.1 ex[i] = 0.0 # Deciding that x is precise! ey[i] = 0.1 # Writing down the error. if (verbose and i < Nverbose) : print " Graph point %2d: x = %6.3f +- %5.3f y = %6.3f +- %5.3f"%(i,x[i], ex[i], y[i], ey[i]) # Define the graph (Npoints, x, y, error x, error y): Graph_x = TGraphErrors(Ndata, x, y, ex, ey) # Plot graph: # ----------- canvas2 = TCanvas( "canvas2", "canvas2", 100, 100, 1200, 600 ) Graph_x.SetTitle("Fit of a graph") Graph_x.GetXaxis().SetRangeUser(-5.0, 5.0) Graph_x.GetXaxis().SetTitle("X") Graph_x.GetYaxis().SetTitle("Y") Graph_x.SetMarkerStyle(20) Graph_x.SetLineColor(kBlue) Graph_x.SetLineWidth(1) Graph_x.Draw("AP") # Fit graph: # ---------- # There are three ways of fitting! # 1: Predefined function (ROOT has gaus, expo, polX, etc.) # 2: Writing function explicitly (for simple functions) # 3: Defining external function (for advanced functions) # Predefined function: fit_x1 = TF1("fit_x1", "pol3", -1.3, 1.3) # "pol3" is the predefined 3rd degree polynomial fit_x1.SetLineColor(kRed) Graph_x.Fit("fit_x1", "R") # Option "R" means "fit in range of function only! # Writing function explicitely: fit_x2 = TF1("fit_x2", "[0] + [1]*x + [2]*x*sin(x) + [3]*x*x*x", -1.3, 1.3) fit_x2.SetParameters(0.0, 2.0, 0.0, 2.0) # Remember to give good (reasonable) starting values! fit_x2.SetLineColor(kMagenta) Graph_x.Fit("fit_x2","R+") # Option "+" adds fit (i.e. not deleting old one!). # Defining external function (which can include "if" statements and everything else!): def func_advanced (x, p) : if (x[0] > 0.0) : return p[0] + p[1]*x[0] + p[2]*x[0]*sin(x[0]) + p[3]*x[0]*x[0]*x[0] + p[4] else : return p[0] + p[1]*x[0] + p[2]*x[0]*sin(x[0]) + p[3]*x[0]*x[0]*x[0] fit_x3 = TF1("fit_x3", func_advanced, -1.3, 1.3, 5) # Here you need to ALSO define number of variables fit_x3.SetParameters(0.0, 2.0, 0.0, 2.0, 0.5) # Remember to give good (reasonable) starting values! fit_x3.SetLineColor(kBlue) Graph_x.Fit("fit_x3", "R+") canvas2.Update() if (SavePlots): canvas2.SaveAs("Graph.pdf") # ---------------------------------------------------------------------------------- raw_input('Press Enter to exit')