ATLAS Offline Software
create_input.py
Go to the documentation of this file.
1 # Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
2 
3 # this file do not work out of the box
4 
5 import shutil
6 import ROOT
7 from array import array
8 import numpy as np
9 
10 try:
11  ROOT.gROOT.ProcessLine(".x ~/atlasstyle/AtlasStyle.C")
12  ROOT.SetAtlasStyle()
13 except AttributeError:
14  pass
15 
16 old_filename = "/home/turra/simon/latest_simon.root"
17 old_file = ROOT.TFile.Open(old_filename)
18 old_scale_histogram = old_file.Get("Scales/es2012c/alphaZee_errStat")
19 old_ct_histogram = old_file.Get("Resolution/es2012c/ctZee_errStat")
20 old_scale_sys = old_file.Get("Scales/es2012c/alphaZee_errSyst")
21 old_ct_sys = old_file.Get("Resolution/es2012c/ctZee_errSyst")
22 
23 assert old_scale_histogram
24 assert old_ct_histogram
25 assert old_scale_sys
26 assert old_ct_sys
27 
28 file_christophe = ROOT.TFile("~/ElectronEnergyScaleFactor.root")
29 scales_christophe = file_christophe.Get("alpha68")
30 ct_christophe = file_christophe.Get("sigma24")
31 # sum of 8 / 13 TeV diff + 64 / 32 bins diff
32 scales_sys_christophe = file_christophe.Get("systAlpha")
33 ct_sys_christophe = file_christophe.Get("systSigma") # 8 / 13 TeV diff
34 
35 assert scales_christophe
36 assert ct_christophe
37 assert scales_sys_christophe
38 assert ct_sys_christophe
39 
40 
41 def qsum_histograms(histo1, histo2):
42  new_histo = histo1.Clone()
43  new_histo.Reset()
44  for ibin in range(histo1.GetNbinsX() + 2):
45  value1 = histo1.GetBinContent(ibin)
46  central_value = histo1.GetBinCenter(ibin)
47  ibin2 = histo2.FindBin(central_value)
48  value2 = histo2.GetBinContent(ibin2)
49  new_value = np.sqrt(value1 * value1 + value2 * value2)
50  new_histo.SetBinContent(ibin, new_value)
51  return new_histo
52 
53 
54 def copy_bin(histo, eta_range):
55  """ avoid bin with value == 0 """
56  empty_bins_in_transition = set()
57  for eta in eta_range:
58  bin = histo.FindBin(eta)
59  if (histo.GetBinContent(bin) == 0):
60  empty_bins_in_transition.add(bin)
61  assert len(empty_bins_in_transition) == 1
62  empty_bin = list(empty_bins_in_transition)[0]
63  if eta_range[-1] > 0:
64  assert(histo_scale.GetBinContent(empty_bin - 1) != 0)
65  histo.SetBinContent(empty_bin, histo.GetBinContent(empty_bin - 1))
66  histo.SetBinError(empty_bin, histo.GetBinError(empty_bin - 1))
67  else:
68  assert(histo_scale.GetBinContent(empty_bin + 1) != 0)
69  histo.SetBinContent(empty_bin, histo.GetBinContent(empty_bin + 1))
70  histo.SetBinError(empty_bin, histo.GetBinError(empty_bin + 1))
71 
72 
73 def copy_dir(source, destination):
74  # equivalent of cp -r source/* destination
75  savdir = ROOT.gDirectory
76  destination.cd()
77  # loop on all entries of this directory
78  for key in source.GetListOfKeys():
79  classname = key.GetClassName()
80  cl = ROOT.gROOT.GetClass(classname)
81  if (not cl):
82  continue
83  if (cl.InheritsFrom(ROOT.TDirectory.Class())):
84  d = destination.mkdir(key.GetName())
85  copy_dir(source.GetDirectory(key.GetName()), d)
86  destination.cd()
87  elif (cl.InheritsFrom(ROOT.TTree.Class())):
88  T = source.Get(key.GetName())
89  destination.cd()
90  newT = T.CloneTree(-1, "fast")
91  newT.Write()
92  elif (cl.InheritsFrom(ROOT.TList.Class())):
93  source.cd()
94  obj = key.ReadObj()
95  destination.cd()
96  obj.Write(obj.GetName(), ROOT.TObject.kSingleKey)
97  del obj
98  else:
99  source.cd()
100  obj = key.ReadObj()
101  destination.cd()
102  obj.Write()
103  del obj
104 
105  destination.SaveSelf(True)
106  savdir.cd()
107 
108 
109 def create_structured_file(path):
110  output_file = ROOT.TFile(path, "recreate")
111  output_file.mkdir("Scales")
112  output_file.mkdir("Resolution")
113  output_file.mkdir("PSRecalibration")
114  output_file.mkdir("S1Recalibration")
115  return output_file
116 
117 
118 def create_new_directories(output_file):
119  output_file.Get("Scales").mkdir("es2015PRE")
120  output_file.Get("Resolution").mkdir("es2015PRE")
121  output_file.Get("PSRecalibration").mkdir("es2015PRE")
122  output_file.Get("S1Recalibration").mkdir("es2015PRE")
123 
124 
125 def merge_histograms(old, new, merge_error=True):
126  new_binning = []
127  new_values = []
128  new_errors = []
129  UNDERFLOW = 0
130  OVERFLOW = new.GetNbinsX() + 1
131 
132  for iold in range(1, old.GetNbinsX()):
133  low = old.GetBinLowEdge(iold)
134  r = low + old.GetBinWidth(iold)
135 
136  il_new = new.FindFixBin(low)
137  ir_new = new.FindFixBin(r)
138  remainer = None
139 
140  if il_new == UNDERFLOW and ir_new == UNDERFLOW:
141  new_binning.append((low, r))
142  new_values.append(old.GetBinContent(iold))
143  new_errors.append(old.GetBinError(iold))
144 
145  elif il_new == UNDERFLOW and ir_new > UNDERFLOW and ir_new < OVERFLOW:
146  if abs(new.GetBinLowEdge(1) - low) < 1E-100:
147  continue
148  new_binning.append((low, new.GetBinLowEdge(1)))
149  new_values.append(old.GetBinContent(iold))
150  new_errors.append(old.GetBinError(iold))
151  if ir_new == OVERFLOW:
152  remainer = iold
153  break
154  last_old = iold
155 
156  for inew in range(1, new.GetNbinsX() + 1):
157  low = new.GetBinLowEdge(inew)
158  r = low + new.GetBinWidth(inew)
159  new_binning.append((low, r))
160  new_values.append(new.GetBinContent(inew))
161  new_errors.append(new.GetBinError(inew))
162 
163  if remainer is not None:
164  new_binning.append((new.GetBinLowEdge(new.GetNbinsX()),
165  old.GetBinLowEdge(remainer)
166  + old.GetBinWidth(remainer)))
167  new_values.append(old.GetBinContent(remainer))
168  new_errors.append(old.GetBinError(remainer))
169 
170  for iold in range(last_old, old.GetNbinsX() + 1):
171  low = old.GetBinLowEdge(iold)
172  r = low + old.GetBinWidth(iold)
173 
174  il_new = new.FindFixBin(low)
175  ir_new = new.FindFixBin(r)
176 
177  if il_new == OVERFLOW and ir_new == OVERFLOW:
178  new_binning.append((low, r))
179  new_values.append(old.GetBinContent(iold))
180  new_errors.append(old.GetBinError(iold))
181  elif il_new < OVERFLOW and ir_new == OVERFLOW:
182  if abs(new.GetBinLowEdge(new.GetNbinsX() + 1) - r) < 1E-100:
183  continue
184  new_binning.append((new.GetBinLowEdge(new.GetNbinsX() + 1), r))
185  new_values.append(old.GetBinContent(iold))
186  new_errors.append(old.GetBinError(iold))
187 
188  new_edges = array('f', [x[0] for x in new_binning] + [new_binning[-1][1]])
189  histo_type = type(new)
190  result = histo_type(new.GetName(), new.GetTitle(),
191  len(new_edges) - 1, new_edges)
192  for i, (v, e) in enumerate(zip(new_values, new_errors), 1):
193  result.SetBinContent(i, v)
194  if merge_error:
195  result.SetBinError(i, e)
196  return result
197 
198 
199 histo_scale = merge_histograms(old_scale_histogram, scales_christophe)
200 histo_ct = merge_histograms(old_ct_histogram, ct_christophe)
201 
202 
203 copy_bin(histo_scale, np.arange(2, 2.7, 0.001))
204 copy_bin(histo_scale, np.arange(-2.7, -2, 0.001))
205 
206 copy_bin(histo_ct, np.arange(2, 2.7, 0.001))
207 copy_bin(histo_ct, np.arange(-2.7, -2, 0.001))
208 
209 
210 histo_scale.SetTitle("Zee energy scale factor - es2015 prerecommendation")
211 histo_ct.SetTitle("Zee additional constant term - es2015 prerecommendation")
212 
213 for h in histo_scale, histo_ct:
214  h.SetMarkerColor(ROOT.kBlack)
215  h.SetMarkerStyle(20)
216 
217 canvas1 = ROOT.TCanvas()
218 histo_scale.Draw("")
219 old_scale_histogram.Draw("sameL")
220 scales_christophe.Draw("sameL")
221 
222 canvas2 = ROOT.TCanvas()
223 histo_ct.Draw("L")
224 old_ct_histogram.Draw("sameL")
225 ct_christophe.Draw("sameL")
226 
227 histo_scale.SetName("alphaZee_errStat")
228 histo_ct.SetName("ctZee_errStat")
229 
230 # this created a file structured as the official one, with empty directories
231 # output_file = create_structured_file("calibration_constants_run2.root")
232 print(old_filename)
233 shutil.copy2(old_filename, "xxx.root")
234 output_file = ROOT.TFile("xxx.root", "update")
235 create_new_directories(output_file)
236 
237 output_file.GetDirectory("Scales").GetDirectory(
238  "es2015PRE").WriteObject(histo_scale, "alphaZee_errStat")
239 output_file.GetDirectory("Resolution").GetDirectory(
240  "es2015PRE").WriteObject(histo_ct, "ctZee_errStat")
241 
242 # systematics
243 
244 new_scale_sys = qsum_histograms(scales_sys_christophe, old_scale_sys)
245 new_ct_sys = qsum_histograms(ct_sys_christophe, old_ct_sys)
246 new_scale_sys = merge_histograms(old_scale_sys, new_scale_sys, False)
247 new_ct_sys = merge_histograms(old_ct_sys, new_ct_sys, False)
248 new_scale_sys.SetTitle(
249  "es2015PRE scales sys = es2012c sys + 7/13 TeV diff + 34/68 bin diff")
250 new_ct_sys.SetTitle("es2015 ct sys = es2012c sys + 7/13 TeV diff")
251 output_file.GetDirectory("Scales").GetDirectory(
252  "es2015PRE").WriteObject(new_scale_sys, "alphaZee_errSyst")
253 output_file.GetDirectory("Resolution").GetDirectory(
254  "es2015PRE").WriteObject(new_ct_sys, "ctZee_errSyst")
255 
256 legend3 = ROOT.TLegend(0.6, 0.7, 0.9, 0.9)
257 legend3.SetBorderSize(0)
258 canvas3 = ROOT.TCanvas()
259 old_scale_sys.GetXaxis().SetTitle("#eta")
260 old_scale_sys.GetYaxis().SetTitle("scale systematics")
261 old_scale_sys.SetLineColor(ROOT.kRed)
262 old_scale_sys.Draw()
263 scales_sys_christophe.SetLineColor(ROOT.kBlue)
264 scales_sys_christophe.Draw("same")
265 new_scale_sys.Draw("same")
266 legend3.AddEntry(scales_sys_christophe, "8/13 TeV diff + 34/68 bin diff", "L")
267 legend3.AddEntry(old_scale_sys, "2012c systematics", "L")
268 legend3.AddEntry(new_scale_sys, "sum (2015PRE)", "L")
269 legend3.Draw()
270 
271 legend4 = ROOT.TLegend(0.6, 0.7, 0.9, 0.9)
272 legend4.SetBorderSize(0)
273 canvas4 = ROOT.TCanvas()
274 old_ct_sys.GetXaxis().SetTitle("#eta")
275 old_ct_sys.GetYaxis().SetTitle("ct systematics")
276 old_ct_sys.SetLineColor(ROOT.kRed)
277 old_ct_sys.Draw()
278 ct_sys_christophe.SetLineColor(ROOT.kBlue)
279 ct_sys_christophe.Draw("same")
280 new_ct_sys.Draw("same")
281 legend4.AddEntry(ct_sys_christophe, "8/13 TeV diff", "L")
282 legend4.AddEntry(old_ct_sys, "2012c systematics", "L")
283 legend4.AddEntry(new_ct_sys, "sum (2015PRE)", "L")
284 legend4.Draw()
285 
286 # stefano input for sensitivity
287 
288 stefano_file = ROOT.TFile("egammaEnergyCorrectionDataMC15.root")
289 copy_dir(stefano_file.GetDirectory("PSRecalibration"),
290  output_file.GetDirectory("PSRecalibration").GetDirectory("es2015PRE"))
291 copy_dir(stefano_file.GetDirectory("S1Recalibration"),
292  output_file.GetDirectory("S1Recalibration").GetDirectory("es2015PRE"))
293 
294 # correction for uA2MeV first week 2015
295 
296 f_ua2mev = ROOT.TFile.Open("~/uA2MeV.root")
297 histo_ua2mev = f_ua2mev.Get("histo_uA2MeV_week12")
298 output_file.Get("Scales").Get("es2015PRE").WriteTObject(histo_ua2mev)
299 
300 input("Press Key")
create_input.copy_bin
def copy_bin(histo, eta_range)
Definition: create_input.py:54
create_input.create_structured_file
def create_structured_file(path)
Definition: create_input.py:109
plot_material.mkdir
def mkdir(path, recursive=True)
Definition: plot_material.py:16
Get
T * Get(TFile &f, const std::string &n, const std::string &dir="", const chainmap_t *chainmap=0, std::vector< std::string > *saved=0)
get a histogram given a path, and an optional initial directory if histogram is not found,...
Definition: comparitor.cxx:179
PlotPulseshapeFromCool.input
input
Definition: PlotPulseshapeFromCool.py:106
plotBeamSpotVxVal.range
range
Definition: plotBeamSpotVxVal.py:195
create_input.qsum_histograms
def qsum_histograms(histo1, histo2)
Definition: create_input.py:41
histSizes.list
def list(name, path='/')
Definition: histSizes.py:38
CxxUtils::set
constexpr std::enable_if_t< is_bitmask_v< E >, E & > set(E &lhs, E rhs)
Convenience function to set bits in a class enum bitmask.
Definition: bitmask.h:232
create_input.create_new_directories
def create_new_directories(output_file)
Definition: create_input.py:118
print
void print(char *figname, TCanvas *c1)
Definition: TRTCalib_StrawStatusPlots.cxx:25
array
python.CaloScaleNoiseConfig.type
type
Definition: CaloScaleNoiseConfig.py:78
create_input.copy_dir
def copy_dir(source, destination)
Definition: create_input.py:73
create_input.merge_histograms
def merge_histograms(old, new, merge_error=True)
Definition: create_input.py:125
Generated on Sun Dec 22 2024 21:08:46 for ATLAS Offline Software by   doxygen 1.8.18