-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathrunSignalFits.py
131 lines (109 loc) · 4.17 KB
/
runSignalFits.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
import ROOT
import sys
import os
import numpy as np
import json
ROOT.gROOT.SetBatch(True)
if __name__ == "__main__":
from modules.utils import parseRuns
run_start, run_end = parseRuns()
from modules.runinfo import runinfo, GetFitRange, IsMuonRun
from modules.fitFunction import runFit, runMIPFit
runs = []
energys = []
mus = []
muEs = []
sigmas = []
sigmaEs = []
mus_linear = []
muEs_linear = []
sigmas_linear = []
sigmaEs_linear = []
for run in range(run_start, run_end+1):
fname = f"calibrated/Run{run}_list.root"
if not os.path.exists(fname):
print(f"File {fname} does not exist")
continue
if run not in runinfo:
print(f"Run {run} not in run info")
continue
if IsMuonRun(run):
print(f"Run {run} is a muon run")
continue
be, hasAtten, hasFilter, _ = runinfo[run]
energy = be * 8.0
print(f"Run {run} has energy {energy} GeV")
print(f"attenuation {hasAtten}, neural density filter {hasFilter}")
fitranges = GetFitRange(int(energy), hasAtten, hasFilter)
if fitranges is None:
print(f"Fit range not found for run {run}")
continue
print(f"Fit ranges: {fitranges}")
f = ROOT.TFile(fname)
hcal = f.Get("hcal_mip")
hcal_linear = f.Get("hcal_linear")
hcal_unc = f.Get("hcal_unc")
if energy == 8.0:
hcal.Rebin(2)
elif energy >= 12.0 and energy < 20.0:
hcal.Rebin(4)
elif energy >= 20.0:
hcal.Rebin(5)
if not hasAtten and not hasFilter:
hcal.Rebin(2)
(mu, muE), (sigma, sigmaE) = runFit(hcal, f"cal_{run}", xmin=fitranges[0], xmax=fitranges[1],
xfitmin=fitranges[2], xfitmax=fitranges[3], outdir="plots/MIPCalibed/Fits/")
# (mu, muE), (sigma, sigmaE) = runMIPFit(hcal, f"cal_{run}", xmin=fitranges[0], xmax=fitranges[1],
# xfitmin=fitranges[2], xfitmax=fitranges[3], outdir="plots/MIPCalibed/MIPFits/")
if hcal_linear:
if energy == 8.0:
hcal_linear.Rebin(2)
elif energy == 12.0:
hcal_linear.Rebin(4)
elif energy >= 20.0:
hcal_linear.Rebin(5)
fit_ranges = GetFitRange(
int(energy), hasAtten, hasFilter, isLinearRegression=True)
(mu_linear, muE_linear), (sigma_linear, sigmaE_linear) = runFit(hcal_linear, f"linear_{run}",
xmin=fitranges[0], xmax=fitranges[1], xfitmin=fitranges[2], xfitmax=fitranges[3], outdir="plots/LinearRegressed/Fits/")
if hcal_unc:
if energy == 8.0:
hcal_unc.Rebin(2)
elif energy == 12.0:
hcal_unc.Rebin(4)
elif energy >= 20.0:
hcal_unc.Rebin(5)
runFit(hcal_unc, f"uncal_{run}", xmin=fitranges[0], xmax=fitranges[1],
xfitmin=fitranges[2], xfitmax=fitranges[3], be=be, hasAtten=hasAtten)
runs.append(run)
energys.append(energy)
mus.append(mu)
muEs.append(muE)
sigmas.append(sigma / mu)
sigmaEs.append(sigmaE / mu)
mus_linear.append(mu_linear)
muEs_linear.append(muE_linear)
sigmas_linear.append(sigma_linear / mu_linear)
sigmaEs_linear.append(sigmaE_linear / mu_linear)
fitresults = {
"runs": runs,
"energys": energys,
"mus": mus,
"muEs": muEs,
"sigmas": sigmas,
"sigmaEs": sigmaEs,
"mus_linear": mus_linear,
"muEs_linear": muEs_linear,
"sigmas_linear": sigmas_linear,
"sigmaEs_linear": sigmaEs_linear
}
if not os.path.exists("results"):
print("Creating results directory")
os.makedirs("results")
with open(f"results/fitresults_Run{run_start}_{run_end}.json", "w") as f:
json.dump(fitresults, f)
print("energys = ", energys)
print("mus = ", mus)
print("muEs = ", muEs)
print("sigmas = ", sigmas)
print("sigmaEs = ", sigmaEs)