Calibration ignores hcal component and tries to fit ecal with linear and fractional term, E^(3/2) to match with HCal/ECal linear extrapolation data. Two calibration constants are gained - a linear and fractional term for the total ecal with assumed relative factor of 2. Done only in a full fit method
Backscatter removed, no tail cuts
In [1]:
using LCIO
using Plots
using Glob
using LsqFit
using StatsBase
using Distributions
using StatPlots
gr();
In [2]:
# read the files for calibration - store files in `singleParticles` directory
fileList = readdir(glob"reco_5000a*0phi*", "singleParticles")
phi=0.0
# build the list of input energies
energyString = r"_(\d+)GeV"
# get (match) the energy string from the filename, convert (parse) from string to Int16, build a dictionary from ints to filenames
energyMap = Dict((parse(Int16, match(energyString, fn)[1]), fn) for fn in fileList)
for (energy, filename) in energyMap
println(energy, " GeV:\t", filename)
end
In [3]:
function getHitsFromFile(filename)
eCalEnergies = Float64[]
hCalEnergies = Float64[]
eCalBackEns = Float64[]
eCalLengths = Int16[]
LCIO.open(filename) do reader
for (idx, event) in enumerate(reader)
if idx > 10000
break
end
eCalBackEn = 0.0
hCalEnergy = 0.0
eCalEnergy = 0.0
eCalLength = 0
# sum up the uncalibrated ECalHits
# this needs to be sorted by layer, so we need a decoder
EcalBarrelHits = getCollection(event, "ECalBarrelHits")
decode = CellIDDecoder(EcalBarrelHits)
for hit in EcalBarrelHits
##eliminate backscatter - if angle > 0.2 rad from center at "phi"
angle=(acos((cos(phi*pi/180.)*getPosition(hit)[1]+sin(phi*pi/180.)*getPosition(hit)[2])/sqrt(getPosition(hit)[1]*getPosition(hit)[1]+getPosition(hit)[2]*getPosition(hit)[2]+getPosition(hit)[3]*getPosition(hit)[3])))
if getPosition(hit)[1]<0
angle=-angle
end
if phi*pi/180+0.2>angle>phi*pi/180-0.2
# calibrate the hits in the later layers with a higher number, because they are behind thicker tungsten slabs
factor = decode(hit)["layer"] < 21 ? 1 : 2
if decode(hit)["layer"]>0
eCalEnergy += factor*getEnergy(hit)
eCalLength += 1
if decode(hit)["layer"]>=29
eCalBackEn+=getEnergy(hit)
end
end
end
end
HcalBarrelHits = getCollection(event,"HCalBarrelHits")
for hhit in HcalBarrelHits
hCalEnergy += getEnergy(hhit) #sum up all raw HCal deposits
end
# fixme: Simple outlier cut
if eCalLength < 30
continue
end
push!(eCalEnergies, eCalEnergy)
push!(eCalLengths, eCalLength)
push!(eCalBackEns, eCalBackEn)
push!(hCalEnergies,hCalEnergy)
end
end
return eCalEnergies, eCalLengths, eCalBackEns, hCalEnergies
end;
In [4]:
eHits = Dict{Int16, Vector{Float64}}()
hHits = Dict{Int16, Vector{Float64}}()
eCount = Dict{Int16, Vector{Int16}}()
eHitsBack = Dict{Int16, Vector{Float64}}()
for (energy, filename) in energyMap
println("Processing file for ", energy, " GeV")
eCal, nEhits, eCalBack, hCal = getHitsFromFile(filename)
eHits[energy] = eCal
hHits[energy] = hCal
eCount[energy] = nEhits
eHitsBack[energy] = eCalBack
end
In [5]:
sigmasUncalib=Vector{Float64}()
meansUncalib=Vector{Float64}()
gaus=Distributions.fit(Normal,eHits[100])
plot(gaus,linewidth=3)
for energy in keys(eHits)
gausFn=Distributions.fit(Normal,eHits[energy]) #normalized by area
println(energy," GeV: μ = ", gausFn.μ, "\tσ = ", gausFn.σ)
push!(sigmasUncalib, gausFn.σ)
push!(meansUncalib, gausFn.μ)
plot!(gausFn,linewidth=3)
end
histogram!([eHits[energy]for energy in keys(eHits)], fillalpha=0.5, linewidth=0, label=map(string, keys(eHits)), xlabel="Uncalibrated Energy [GeV]",normalize=true)
xaxis!(:log10)
Out[5]:
In [6]:
# this function attempts a global fit and minimizes the offset b of the quadratic form y=axx+mx+b
# parameters are ecal energies (×2 for the hits in the outer layers), hcal energies, particle energy
function lineFitter(ecal, truValues)
function model(x, p)
energies = Dict{Int16, Float64}()
for e in truValues
calibrated=p[1].*x[e] + p[2].*x[e].^(3.0/2.0)
n=Distributions.fit(Normal,calibrated)
energies[e] = n.μ
end
# we are optimizing for the ratio of reconstructed energies to true values
return [energies[e]/e for e in truValues]
end
fit = curve_fit(model, ecal, 1.0, [50.0,0.0])
errors=estimate_errors(fit,0.95)
return fit.param,errors
end
ECal, ECalErr = lineFitter(eHits, keys(eHits))
println("ECal calibration constant: ", ECal[1], " +/- ", ECalErr[1])
println("ECal fractional calibration constant: ", ECal[2], " +/- ", ECalErr[2])
In [7]:
sigmasCalib = Vector{Float64}()
meansCalib = Vector{Float64}()
gaus=Distributions.fit(Normal,ECal[2] .* eHits[100].^(3.0/2.0) + ECal[1] .* eHits[100])
plot(gaus,linewidth=3,label="")
for energy in keys(eHits)
gausFn=Distributions.fit(Normal,ECal[2] .* eHits[energy].^ (3.0/2.0) + ECal[1] .* eHits[energy])
println(energy," GeV: μ = ", gausFn.μ, "\tσ = ", gausFn.σ)
push!(sigmasCalib,gausFn.σ)
push!(meansCalib,gausFn.μ)
plot!(gausFn,linewidth=3,label="")
end
histogram!([ECal[2] .* eHits[energy].^(3.0/2.0) + ECal[1] .* eHits[energy] for energy in keys(eHits)], fillalpha=0.5, linewidth=0, label=map(string, keys(eHits)),normalize=true,xlabel="Calibrated Energy [GeV]")
xaxis!(:log10)
Out[7]:
In [8]:
calibrated = Vector{Float64}()
x = Vector{Float64}()
xLin=[1.,2.,5.,10.,20.,50.,100.]
for e in keys(eHits)
d = ECal[1] .* eHits[e] + ECal[2] .* eHits[e] .^ (3.0/2.0)
gauss=Distributions.fit(Normal,d)
push!(calibrated, gauss.μ)
push!(x, 1.0*e)
end
line(x, p) = p[1] + p[2]*x
l = curve_fit(line, x, calibrated, [0.0, 1.0])
plot(x, calibrated, yerr=sigmasCalib, marker=stroke(2), line=false, label="data",xlabel="Initial Energy [GeV]",ylabel="Calibrated Energy [GeV]")
leg = @sprintf("y=%.2f + %.2f*x ", l.param[1], l.param[2])
plot!(x, l.param[1]+l.param[2]*x, label=leg)
plot!(xLin,xLin,label="Diagonal")
xaxis!(:log10,(0.9,110.0))
yaxis!(:log10,(0.9,110.0))
Out[8]:
In [9]:
calibratedLin = Vector{Float64}()
for e in keys(eHits)
d = ECal[1] .* eHits[e]
gauss=Distributions.fit(Normal,d)
push!(calibratedLin, gauss.μ)
end
l2 = curve_fit(line, x, calibratedLin, [0.0, 1.0])
plot(x, calibrated, yerr=sigmasCalib, marker=stroke(2), line=false, label="data",xlabel="Initial Energy [GeV]",ylabel="Linearly Calibrated Energy [GeV]")
leg2 = @sprintf("y=%.2f + %.2f*x ", l2.param[1], l2.param[2])
plot!(x, l2.param[1]+l2.param[2]*x, label=leg2)
plot!(xLin,xLin,label="Diagonal")
xaxis!(:log10,(0.9,110.0))
yaxis!(:log10,(0.9,110.0))
Out[9]:
In [10]:
linCalibRatio = Vector{Float64}()
linCalibRatioErr = Vector{Float64}()
for e in keys(eHits)
d = (eHits[e]*ECal[1]) / (1.0*e)
gauss=Distributions.fit(Normal,d)
push!(linCalibRatio, gauss.μ)
push!(linCalibRatioErr, gauss.σ/sqrt(5000))
end
plot(x, linCalibRatio, yerr=linCalibRatioErr, marker=stroke(2), markersize=10, markercolor=:blue,line=false, label="Contained ",xlabel="Initial Energy [GeV]",ylabel="Energy Term / Initial Energy")
xaxis!(:log10,(0.9,110.0))
Out[10]:
In [11]:
calibRatio = Vector{Float64}()
calibRatioErr = Vector{Float64}()
for e in keys(eHits)
d = (ECal[1] .* eHits[e] +ECal[2] .* eHits[e] .^(3.0/2.0))./ (1.0.*e)
gauss=Distributions.fit(Normal,d)
push!(calibRatio, gauss.μ)
push!(calibRatioErr, gauss.σ/sqrt(5000))
end
plot(x, linCalibRatio, yerr=linCalibRatioErr, marker=stroke(2), markersize=10, markercolor=:blue,line=false, label="Contained")
plot!(x, calibRatio, yerr=calibRatioErr, marker=stroke(2),markersize=10, markercolor=:red,line=false,label="Corrected For Leakage Based on Shower Shape ",xlabel="Initial Energy [GeV]",ylabel="Energy Term / Initial Energy")
xaxis!(:log10,(0.9,110.0))
Out[11]:
In [12]:
enRes = Vector{Float64}()
for index in 1:7
push!(enRes,sigmasCalib[index]./meansCalib[index])
end
plot(x, abs((1.0-calibRatio)./enRes), yerr=(1-calibRatio)./enRes./sqrt(5000),marker=stroke(2), markersize=10, markercolor=:salmon,line=false, leg=false,xlabel="Initial Energy [GeV]",ylabel="(1-Calibrated Energy / Initial Energy)/Energy Resolution")
xaxis!(:log10,(0.9,110.0))
yaxis!((-0.001,0.03))
Out[12]:
In [13]:
plot(x, abs((1.0-linCalibRatio)./enRes), yerr=(1-linCalibRatio)./enRes./sqrt(5000),marker=stroke(2),markersize=10, line=false, label="Linear Calibration",xlabel="Initial Energy [GeV]",ylabel="(1-Energy / Initial Energy) / Energy Resolution")
plot!(x, abs((1.0-calibRatio)./enRes), yerr=(1-calibRatio)./enRes./sqrt(5000),marker=stroke(2), markersize=10,label="Full Nonlinear Calibration ", line=false)
xaxis!(:log10,(0.9,110.0))
yaxis!((-0.02,0.65))
Out[13]:
In [14]:
nonlinCalibRatio = Vector{Float64}()
nonlinCalibRatioErr = Vector{Float64}()
for e in keys(eHits)
d = (ECal[2] .* eHits[e] .^(3.0/2.0))./ (1.0.*e)
gauss=Distributions.fit(Normal,d)
push!(nonlinCalibRatio, gauss.μ)
push!(nonlinCalibRatioErr, gauss.σ./sqrt(5000.0))
end
plot(x, linCalibRatio, yerr=linCalibRatioErr, marker=stroke(2), markersize=10, markercolor=:yellow, line=false, label="Linear Calibration Term ",xlabel="Initial Energy [GeV]",ylabel="Energy Term / Initial Energy")
plot!(x, nonlinCalibRatio, yerr=nonlinCalibRatioErr, marker=stroke(2), markersize=10, markercolor=:magenta, line=false,label="1/2 Calibration Term ")
xaxis!(:log10,(0.9,110.0))
yaxis!(:log10,(0.0001,12))
Out[14]:
In [15]:
histogram([eHitsBack[energy]./eHits[energy] for energy in [1,100,10]], bins=(0.0:0.001:0.03),fillalpha=0.5, linewidth=0, label=map(string, keys(eHits)),xlabel="Uncalibrated Deposits in Last 2 Layers [%]")
xaxis!((0.,0.022))
yaxis!((0.,5010.))
Out[15]:
In [16]:
back = Vector{Float64}()
for e in keys(eHits)
push!(back, mean(eHitsBack[e],1)[1]./e.*100)
end
plot(x, nonlinCalibRatio, yerr=nonlinCalibRatioErr, marker=stroke(2),markersize=10, markercolor=:magenta, line=false,label="1/2 Term ",xlabel="Initial Energy [GeV]",right_margin=30px)
plot!(x,back,marker=stroke(2),markersize=10, markercolor=:lightGreen,line=false,label="Uncalibrated Energy in Last 2 ECal Layers / Initial Energy*100 ",ymirror=true,xlink=2, inset=(1,bbox(0,0,1,1)), bg_inside=RGBA(0,0,0,0))
xaxis!(:log10,(0.9,110.0),true)
yaxis!(:log10,(0.001,0.03),true)
Out[16]:
In [17]:
hcalHitsRatio = Vector{Float64}()
for e in keys(eHits)
push!(hcalHitsRatio, mean(hHits[e],1)[1]./e.*50.)
end
plot(x, nonlinCalibRatio, yerr=nonlinCalibRatioErr, marker=stroke(2), markersize=10, markercolor=:magenta,line=false,label="1/2 Term ",xlabel="Initial Energy [GeV]",right_margin=30px)
plot!(x,hcalHitsRatio,marker=stroke(2),markersize=10, markercolor=:orange,line=false,label="Uncalibrated HCal Hits/Initial Energy*50 ",ymirror=true,xlink=2, inset=(1,bbox(0,0,1,1)), bg_inside=RGBA(0,0,0,0))
xaxis!(:log10,(0.9,110.0))
yaxis!(:log10,(0.001,0.03))
Out[17]:
In [18]:
#from MIP vs scaled radiation length plot, linear fit extrapolation
fitIntegralL=[0.0009495162208586111, 0.001249982299277269, 0.0016381364639927664, 0.0023597250460691034, 0.0030983044984638354, 0.004435031937199862, 0.006134404577683722]
#from MIP vs radiation length plot, gamma fit extrapolation
fitIntegralG=[0.00021193102439427688, 0.000340690533861749, 0.0007747472067739187, 0.0009963353823439938, 0.0015979182805538118, 0.0024471932434884614, 0.003817725540069536]
plot(x, nonlinCalibRatio, yerr=nonlinCalibRatioErr, marker=stroke(2), markersize=10, markercolor=:magenta,line=false,label="1/2 Term ",xlabel="Initial Energy [GeV]",right_margin=30px)
plot!(xLin,fitIntegralG,marker=stroke(2),markersize=10, markercolor=:green,line=false,label="Integral value of Gamma Fit Leakage Extrapolation ",ymirror=true,xlink=2, inset=(1,bbox(0,0,1,1)), bg_inside=RGBA(0,0,0,0))
plot!(xLin,fitIntegralL,marker=stroke(2),markersize=10, markercolor=:lightBlue,line=false,label="Integral value of Linear Fit Leakage Extrapolation")
xaxis!(:log10,(0.9,110.0))
yaxis!(:log10,(0.0001,0.05))
Out[18]:
In [19]:
plot(x, nonlinCalibRatio, yerr=nonlinCalibRatioErr, marker=stroke(2),markersize=10, markercolor=:magenta, line=false,label="1/2 Term ",xlabel="Initial Energy [GeV]",right_margin=30px)
plot!(x,back,marker=stroke(2),markersize=10, markercolor=:lightGreen,line=false,label="Uncalibrated Energy in Last 2 Layers / Initial Energy*100 ")
plot!(x,hcalHitsRatio,marker=stroke(2),markersize=10, markercolor=:orange,line=false,label="Uncalibrated HCal Hits/Initial Energy*50 ")
plot!(xLin,fitIntegralG,marker=stroke(2),markersize=10, markercolor=:green,line=false,label="Integral value of Gamma Fit Leakage Extrapolation")
plot!(xLin,fitIntegralL,marker=stroke(2),markersize=10, markercolor=:lightBlue,line=false,label="Integral value of Linear Fit Leakage Extrapolation ",ymirror=true,xlink=2, inset=(1,bbox(0,0,1,1)), bg_inside=RGBA(0,0,0,0))
xaxis!(:log10,(0.9,110.0))
yaxis!(:log10,(0.00008,0.1))
Out[19]:
In [20]:
hcalHitsScaled = Vector{Float64}()
nonlinCalibRatioScaled = Vector{Float64}()
nonlinCalibRatioScaledErr = Vector{Float64}()
integralLScaled = Vector{Float64}()
integralGScaled = Vector{Float64}()
backScaled = Vector{Float64}()
for e in keys(eHits)
d = (ECal[2] .* eHits[e] .^(3.0/2.0))./ (1.0.*e) ./ sqrt(1.0.*e)
gauss=Distributions.fit(Normal,d)
push!(nonlinCalibRatioScaled, gauss.μ)
push!(nonlinCalibRatioScaledErr, gauss.σ./sqrt(5000.0))
end
for e in keys(eHits)
push!(hcalHitsScaled, mean(hHits[e],1)[1]./e/sqrt(1.0.*e))
end
for index in 1:7
push!(integralGScaled, fitIntegralG[index]./sqrt(xLin[index]))
push!(integralLScaled, fitIntegralL[index]./sqrt(xLin[index]))
push!(backScaled, back[index]./sqrt(x[index]))
end
plot(x, nonlinCalibRatioScaled, yerr=nonlinCalibRatioScaledErr, marker=stroke(2), markersize=10, markercolor=:darkMagenta,line=false,label="1/2 Term/sqrt(E) ")
plot!(x,hcalHitsScaled,marker=stroke(2),markersize=10, markercolor=:brown,line=false,label="Uncalibrated HCal Hits/Initial Energy/sqrt(E) ")
plot!(x,backScaled,marker=stroke(2),markersize=10, markercolor=:darkCyan,line=false,label="Uncalibrated Energy in Last 2 Layers / Initial Energy *100 /sqrt(E) ")
plot!(xLin,integralGScaled,marker=stroke(2),markersize=10, markercolor=:darkGreen,line=false,label="Integral value of Gamma Fit Leakage Extrapolation/sqrt(E) ")
plot!(xLin,integralLScaled,marker=stroke(2),markersize=10, markercolor=:darkBlue,line=false,label="Integral value of Linear Fit Leakage Extrapolation/sqrt(E) ",xlabel="Initial Energy [GeV]",right_margin=30px,ymirror=true,xlink=2, inset=(1,bbox(0,0,1,1)), bg_inside=RGBA(0,0,0,0))
xaxis!(:log10,(0.9,110.0))
yaxis!(:log10,(0.00001,0.05))
Out[20]:
In [21]:
plot(xLin,fitIntegralG,marker=stroke(2),markersize=10, markercolor=:green,line=false,label="Integral value of Gamma Fit Leakage Extrapolation ")
plot!(xLin,fitIntegralL,marker=stroke(2),markersize=10, markercolor=:lightBlue,line=false,label="Integral value of Linear Fit Leakage Extrapolation ",xlabel="Initial Energy [GeV]",right_margin=30px,ymirror=true,xlink=2, inset=(1,bbox(0,0,1,1)), bg_inside=RGBA(0,0,0,0))
plot!(x,hcalHitsRatio,marker=stroke(2),markersize=10, markercolor=:orange,line=false,label="Uncalibrated HCal Hits/Initial Energy*50 ")
xaxis!(:log10,(0.9,110.0))
yaxis!(:log10,(0.0001,0.05))
Out[21]:
In [ ]: