#!/usr/bin/python
#
# harptool.cgi - cgi script for running the Hall D coherent bremsstrahlung
# harp scan fitter, a ROOT tool written by R.T.Jones.
#
import sys
sys.path.insert(0, "/usr/local/root/lib/root")
import numpy
from ROOT import *
gROOT.IsBatch()
import os
import cgi
import cgitb
cgitb.enable()
slimits = (100, 200)
def print_head():
print "Content-Type: text/html"
print
print ""
print ""
print "Hall D Electron Beam Harp Scan Fitting Tool "
print ""
print ""
print ""
print ""
print ""
def set_parameter(par, desc, default, unit):
if par in form:
value = form.getfirst(par)
else:
value = default
if str("default " + par) in form:
value = default
print "" + desc + ": ",
print " "
def fit_model(var, par):
"""
var[0] = s coordinate (m)
par[0] = s coordinate (m) of focus
par[1] = sigma at the focus (mm)
par[2] = emittance (mm.mrad)
"""
return (par[1]**2 + (par[2]/par[1] * (var[0] - par[0]))**2)**0.5
def fit_model_gradient(var, par):
"""
var[0] = s coordinate (m)
par[0] = s coordinate (m) of focus
par[1] = sigma at the focus (mm)
par[2] = emittance (mm.mrad)
"""
y = (par[1]**2 + (par[2]/par[1] * (var[0] - par[0]))**2)**0.5
dyds0 = -(par[2]/par[1])**2 * (var[0] - par[0]) / y
dydsig0 = (par[1]**2 - (par[2]/par[1] * (var[0] - par[0]))**2) / (y * par[1])
return (dyds0, dydsig0)
def fit_and_plot(sx, sigx, sigxerr, sy, sigy, sigyerr):
zero = [0] * len(sx)
xfitf = TF1("xfitf", fit_model, slimits[0], slimits[1], 3)
xfitf.SetParameter(0, slimits[1])
xfitf.SetParameter(1, 1.0)
xfitf.FixParameter(2, float(form.getfirst("emittance_x")))
xfitf.SetLineColor(kRed)
xdata = TGraphErrors(len(sx), numpy.array(sx, dtype=float),
numpy.array(sigx, dtype=float),
numpy.array(zero, dtype=float),
numpy.array(sigxerr, dtype=float))
xfit = xdata.Fit(xfitf, "qs")
xdata.SetLineColor(kRed)
xdata.SetLineWidth(2)
xdata.SetMarkerColor(kRed)
xdata.SetMarkerStyle(20)
ss = []
x0 = []
x1 = []
par = [xfit.Parameter(0), xfit.Parameter(1), xfit.Parameter(2)]
covar = xfit.GetCovarianceMatrix()
minvar = 0.01**2
for n in range(0, 101):
s = slimits[0] + n * (slimits[1] - slimits[0]) / 100
xmu = fit_model([s], par)
gra = fit_model_gradient([s], par)
ss.append(s)
xsigma = (covar(0,0) * gra[0]**2 + covar(1,1) * gra[1]**2 +
2 * covar(0,1) * gra[0] * gra[1] + minvar)**0.5
x0.append(xmu - xsigma)
x1.append(xmu + xsigma)
if x0[-1] < 0:
x0[-1] = 0
xshade = TGraph(2 * len(ss), numpy.array(ss + ss[::-1], dtype=float),
numpy.array(x0 + x1[::-1], dtype=float))
xshade.SetFillColorAlpha(kRed, 0.2);
yfitf = TF1("yfitf", fit_model, slimits[0], slimits[1], 3)
yfitf.SetParameter(0, slimits[1])
yfitf.SetParameter(1, 1.0)
yfitf.FixParameter(2, float(form.getfirst("emittance_y")))
yfitf.SetLineColor(kBlue)
ydata = TGraphErrors(len(sy), numpy.array(sy, dtype=float),
numpy.array(sigy, dtype=float),
numpy.array(zero, dtype=float),
numpy.array(sigyerr, dtype=float))
yfit = ydata.Fit(yfitf, "qs")
ydata.SetLineColor(kBlue)
ydata.SetLineWidth(2)
ydata.SetMarkerColor(kBlue)
ydata.SetMarkerStyle(20)
ss = []
y0 = []
y1 = []
par = [yfit.Parameter(0), yfit.Parameter(1), yfit.Parameter(2)]
covar = yfit.GetCovarianceMatrix()
minvar = 0.01**2
for n in range(0, 101):
s = slimits[0] + n * (slimits[1] - slimits[0]) / 100
ymu = fit_model([s], par)
gra = fit_model_gradient([s], par)
ss.append(s)
ysigma = (covar(0,0) * gra[0]**2 + covar(1,1) * gra[1]**2 +
2 * covar(0,1) * gra[0] * gra[1] + minvar)**0.5
y0.append(ymu - ysigma)
y1.append(ymu + ysigma)
if y0[-1] < 0:
y0[-1] = 0
yshade = TGraph(2 * len(ss), numpy.array(ss + ss[::-1], dtype=float),
numpy.array(y0 + y1[::-1], dtype=float))
yshade.SetFillColorAlpha(kBlue, 0.2);
c1 = TCanvas("c1", "", 800, 600)
xfitf.SetTitle("sigma x vs accelerator s")
xfitf.GetXaxis().SetTitle("accelerator s coordinate (m)")
xfitf.GetYaxis().SetTitle("#sigma_{x} (mm)")
xfitf.SetMinimum(0)
xfitf.SetMaximum(max(x1))
xfitf.Draw()
xdata.Draw("P")
xshade.Draw("f")
scol = float(form.getfirst("collimator_spos"))
gcol = TGraph(2, numpy.array([scol] * 2, dtype=float), numpy.array([0, 1.5], dtype=float))
gcol.Draw("L")
c1.Update()
c1.Print(workdir + "harp-x-" + fitimage)
yfitf.SetTitle("sigma y vs accelerator s")
yfitf.GetXaxis().SetTitle("accelerator s coordinate (m)")
yfitf.GetYaxis().SetTitle("#sigma_{y} (mm)")
yfitf.SetMinimum(0)
yfitf.SetMaximum(max(y1))
yfitf.Draw()
ydata.Draw("P")
yshade.Draw("f")
gcol.Draw("L")
c1.Update()
c1.Print(workdir + "harp-y-" + fitimage)
print ""
print ""
print ""
print " "
set_parameter("harp5C11_xsigma", "5C11 harp x sigma", 0, "mm")
print " "
print ""
set_parameter("harp5C11_ysigma", "5C11 harp y sigma", 0, "mm")
print " "
print ""
set_parameter("harp5C11B_xsigma", "5C11B harp x sigma", 0, "mm")
print " "
print ""
set_parameter("harp5C11B_ysigma", "5C11B harp y sigma", 0, "mm")
print " "
print ""
set_parameter("radHarp_xsigma", "radiator harp x sigma", 0, "mm")
print " "
print ""
set_parameter("radHarp_ysigma", "radiator harp y sigma", 0, "mm")
print " "
print ""
print ""
print "Enter measured values in the above input boxes, "
print " "
print ""
print ""
print "defaults are generally ok for the other fields. "
print " "
print ""
set_parameter("emittance_x", "x emittance of e-beam", 0.0062, "mm.mrad")
print " "
print ""
set_parameter("emittance_y", "y emittance of e-beam", 0.0044, "mm.mrad")
print " "
print ""
print " "
print ""
print "Insufficient data, please add at least 2 valid measurements with errors and try again!"
print " "
print ""
print "Invalid data, please correct errors in the above data and try again!"
print "