ATLAS Offline Software
zlumi_alleff.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 # Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3 
4 import math
5 
6 def trig_tag_and_probe(h, lb_minus_one_trig_eff):
7  yld = (h[2], h[3])
8  ylderr = (h.GetBinError(2), h.GetBinError(3))
9  A, B = yld
10 
11  if B == 0:
12  return lb_minus_one_trig_eff[0], lb_minus_one_trig_eff[1]
13 
14  eff = 1./(float(A)/B/2.+1.)
15  inverrsq = ((1/2./B)*ylderr[0])**2+((A/2./B**2)*ylderr[1])**2
16  return eff, (inverrsq**.5)*(eff**2)
17 
18 
19 def reco_tag_and_probe(hmo, hms, hno, hns, lb_minus_one_reco_eff, MC = False):
20  bin1 = hmo.GetXaxis().FindBin(86000)
21  bin2 = hmo.GetXaxis().FindBin(95000)
22  matchos, matchoserr = extract(hmo, bin1, bin2)
23  matchss, matchsserr = extract(hms, bin1, bin2)
24  nomatchos, nomatchoserr = extract(hno, bin1, bin2)
25  nomatchss, nomatchsserr = extract(hns, bin1, bin2)
26 
27  if MC:
28 
29  matchss = 0
30  nomatchss = 0
31  matchsserr = 0
32  nomatchsserr = 0
33 
34  A = float(matchos-matchss)
35  Aerr = (matchoserr**2+matchsserr**2)**.5
36  B = float(nomatchos-nomatchss)
37  Berr = (nomatchoserr**2+nomatchsserr**2)**.5
38 
39  if Berr == 0:
40  Berr = 1.
41  if A == 0 or B/A == -1:
42  eff = lb_minus_one_reco_eff[0]
43  err = lb_minus_one_reco_eff[1]
44  else:
45  eff = 1./(1.+B/A)
46  inverrsq = ((-B/A**2)*Aerr)**2+((1./A)*Berr)**2
47  err = (inverrsq**.5)*(eff**2)
48 
49  return eff, err
50 
51 
52 def template_method(hmo, hms, hno, hns, hto, hts, do_scale, lb_minus_one_reco_eff, corr = None, MC = False):
53  no_sign = False
54 
55  tbin1 = hmo.GetXaxis().FindBin(75000)
56  tbin2 = hmo.GetXaxis().FindBin(104000)
57  tbin3 = hmo.GetXaxis().FindBin(120000)
58  tbin4 = hmo.GetXaxis().FindBin(250000)
59 
60  if no_sign:
61  hmo.Add(hms)
62  hno.Add(hns)
63  hto.Add(hts)
64 
65  matchos_peak, matchos_peakerr = extract(hmo, tbin1, tbin2)
66  matchos_tail, matchos_tailerr = extract(hmo, tbin3, tbin4)
67  matchss_peak, matchss_peakerr = extract(hms, tbin1, tbin2)
68  matchss_tail, matchss_tailerr = extract(hms, tbin3, tbin4)
69 
70  nomatchos_peak, nomatchos_peakerr = extract(hno, tbin1, tbin2)
71  nomatchos_tail, nomatchos_tailerr = extract(hno, tbin3, tbin4)
72  nomatchss_peak, nomatchss_peakerr = extract(hns, tbin1, tbin2)
73  nomatchss_tail, nomatchss_tailerr = extract(hns, tbin3, tbin4)
74 
75  templateos_peak, templateos_peakerr = extract(hto, tbin1, tbin2)
76  templateos_tail, templateos_tailerr = extract(hto, tbin3, tbin4)
77  templatess_peak, templatess_peakerr = extract(hts, tbin1, tbin2)
78  templatess_tail, templatess_tailerr = extract(hts, tbin3, tbin4)
79 
80  totalos_peak = matchos_peak + nomatchos_peak
81  totalos_tail = matchos_tail + nomatchos_tail
82  # totalss_peak = matchss_peak + nomatchss_peak
83  # totalss_tail = matchss_tail + nomatchss_tail
84 
85  # totalos_peakerr = pow(totalos_peak, 0.5)
86  # totalss_peakerr = pow(totalss_peak, 0.5)
87  # totalos_tailerr = pow(totalos_tail, 0.5)
88  # totalss_tailerr = pow(totalss_tail, 0.5)
89 
90  if MC:
91  n1 = matchos_peak
92  d1 = totalos_peak
93  sum_w_m, sum_w2_m = extract(hmo, tbin1, tbin2)
94  sum_w_n, sum_w2_n = extract(hno, tbin1, tbin2)
95 
96  sum_w = sum_w_m + sum_w_n
97  sum_w2 = (sum_w2_m**2) + (sum_w2_n**2)
98 
99  eff = n1/d1
100 
101  err = (eff*(1-eff)*sum_w2/sum_w**2)**0.5
102  return eff, err
103 
104  if templatess_tail == 0:
105  eff = lb_minus_one_reco_eff[0]
106  err = lb_minus_one_reco_eff[1]
107  elif templateos_tail == 0:
108  eff = lb_minus_one_reco_eff[0]
109  err = lb_minus_one_reco_eff[1]
110 
111  if templatess_tail != 0 and templateos_tail != 0:
112  n1 = matchos_peak
113  n2 = templateos_peak*(matchss_tail/templatess_tail)
114  d1 = totalos_peak
115  if not do_scale:
116  d2 = templateos_peak*(nomatchos_tail/templateos_tail)
117  else:
118  d2 = templateos_peak*( (totalos_tail - matchos_tail/corr) /templateos_tail )
119 
120  try:
121  eff = (n1 - n2)/(d1 - d2)
122  err = template_method_error(hmo, hms, hno, hns, hto, hts)
123  except ZeroDivisionError:
124  eff = lb_minus_one_reco_eff[0]
125  err = lb_minus_one_reco_eff[1]
126 
127  return eff, err
128 
129 
130 def template_method_error(hmo, hms, hno, hns, hto, hts):
131  bin1 = hmo.GetXaxis().FindBin(75000)
132  bin2 = hmo.GetXaxis().FindBin(104000)
133  bin3 = hmo.GetXaxis().FindBin(120000)
134  bin4 = hmo.GetXaxis().FindBin(250000)
135 
136  a, da = extract(hmo, bin1, bin2)
137  b, db = extract(hms, bin3, bin4)
138  c, dc = extract(hto, bin1, bin2)
139  d, dd = extract(hts, bin3, bin4)
140  e, de = extract(hno, bin1, bin2)
141  f, df = extract(hno, bin3, bin4)
142  g, dg = extract(hto, bin3, bin4)
143 
144  try:
145  dva = 1/(a-(c*f)/g+e)-(a-(b*c)/d)/(a-(c*f)/g+e)**2
146  except ZeroDivisionError:
147  return 1
148  try:
149  dvb = -c/(d*(-(c*f)/g+a+e))
150  except ZeroDivisionError:
151  return 1
152  try:
153  dvc = -(g*((a+e)*b*g-a*d*f))/(d*(f*c+(-a-e)*g)**2)
154  except ZeroDivisionError:
155  return 1
156  try:
157  dvd = (b*c)/((-(c*f)/g+a+e)*d**2)
158  except ZeroDivisionError:
159  return 1
160  try:
161  dve = -(a-(b*c)/d)/(e-(c*f)/g+a)**2
162  except ZeroDivisionError:
163  return 1
164  try:
165  dvf = (c*(a-(b*c)/d))/(g*(-(c*f)/g+a+e)**2)
166  except ZeroDivisionError:
167  return 1
168  try:
169  dvg = -(c*(a*d-b*c)*f)/(d*((a+e)*g-c*f)**2)
170  except ZeroDivisionError:
171  return 1
172 
173  verbose = False
174  if verbose:
175  print(a, da, dva, da*dva)
176  print(b, db, dvb, db*dvb)
177  print(c, dc, dvc, dc*dvc)
178  print(d, dd, dvd, dd*dvd)
179  print(e, de, dve, de*dve)
180  print(f, df, dvf, df*dvf)
181  print(g, dg, dvg, dg*dvg)
182 
183  err = math.sqrt((da*dva)**2 + (db*dvb)**2 + (dc*dvc)**2 + (dd*dvd)**2 + (de*dve)**2 + (df*dvf)**2 + (dg*dvg)**2)
184  return err
185 
186 def container_efficiency(h_photon, h_photon_total_input, h_fit, h_pass, h_tos, h_tss):
187  h_photon_total = h_photon_total_input.Clone()
188  h_temp = h_tos.Clone()
189  h_temp.Add(h_tss)
190 
191  h_bkg = h_fit.Clone()
192  h_bkg.Reset()
193  for xbin in range(1, h_bkg.GetNbinsX()+1):
194  mass = h_bkg.GetBinCenter(xbin)
195  h_bkg.SetBinContent(xbin, h_fit.GetFunction("pol2").Eval(mass))
196  # This is not correct, but use for now as a proxy -> not sure how to return error from TF1
197  if h_fit.GetFunction("pol2").Eval(mass) < 0:
198  print("Negative fit!", xbin, h_fit.GetFunction("pol2").Eval(mass))
199  h_bkg.SetBinError(xbin, 0)
200  else:
201  h_bkg.SetBinError(xbin, pow(h_fit.GetFunction("pol2").Eval(mass), 0.5))
202 
203 
204  bkg, bkgerr = extract(h_bkg, h_bkg.FindBin(86000), h_bkg.FindBin(96000))
205 
206  bin1 = h_pass.FindBin(120000)
207  bin2 = h_pass.FindBin(250000)
208 
209  bin86 = h_pass.FindBin(86000)
210  bin96 = h_pass.FindBin(96000)
211 
212  xmin = 1
213  xmax = h_photon.GetNbinsX()
214 
215  a, da = extract(h_pass, bin86, bin96)
216  b, db = extract(h_temp, bin86, bin96)
217  c, dc = extract(h_photon_total, bin86, bin96)
218  d, dd = bkg, bkgerr
219 
220  s1, ds1 = extract(h_pass, bin1, bin2)
221  s2, ds2 = extract(h_temp, bin1, bin2)
222  s3, ds3 = extract(h_photon, xmin, xmax)
223  s4, ds4 = extract(h_photon_total, xmin, xmax)
224 
225  numer = a - (b * s1/s2)
226  denom = s3/s4 * (c - d)
227 
228  v2 = pow(numer + denom, 2)
229  try:
230  dda = denom/v2
231  except ZeroDivisionError:
232  dda = 1
233  try:
234  ddb = (s1/s2)*(-denom)/v2
235  except ZeroDivisionError:
236  ddb = 1
237  try:
238  ddc = -(s3/s4) * numer/v2
239  except ZeroDivisionError:
240  ddc = 1
241  try:
242  ddd = (s3/s4) * numer/v2
243  except ZeroDivisionError:
244  ddd = 1
245  try:
246  dds1 = (denom * b/s2)/v2
247  except ZeroDivisionError:
248  dds1 = 1
249  try:
250  dds2 = (denom * b * s1/s2**2)/v2
251  except ZeroDivisionError:
252  dds2 = 1
253  try:
254  dds3 = -(numer) * (denom)/s3 /v2
255  except ZeroDivisionError:
256  dds3 = 1
257  try:
258  dds4 = -(numer) * (denom)/s4 /v2
259  except ZeroDivisionError:
260  dds4 = 1
261 
262  try:
263  conteff = numer / (numer + denom)
264  conteffstat = math.sqrt((da*dda)**2 + (db*ddb)**2 + (dc*ddc)**2 + (dd*ddd)**2 + (ds1*dds1)**2 + (ds2*dds2)**2 + (ds3*dds3)**2 + (ds4*dds4)**2)
265  except Exception:
266  conteff = 0.0
267  conteffstat = 0.0
268 
269  return conteff, conteffstat
270 
271 
272 def extract(histogram, bin1, bin2):
273  error = 0.0
274  for xbin in range(bin1, bin2):
275  error += pow(histogram.GetBinError(xbin), 2)
276  error = math.sqrt(error)
277 
278  return (histogram.Integral(bin1, bin2), error)
279 
tools.zlumi_alleff.container_efficiency
def container_efficiency(h_photon, h_photon_total_input, h_fit, h_pass, h_tos, h_tss)
Definition: zlumi_alleff.py:186
RootHelpers::FindBin
Int_t FindBin(const TAxis *axis, const double x)
Definition: RootHelpers.cxx:14
conifer::pow
constexpr int pow(int x)
Definition: conifer.h:20
tools.zlumi_alleff.template_method
def template_method(hmo, hms, hno, hns, hto, hts, do_scale, lb_minus_one_reco_eff, corr=None, MC=False)
Definition: zlumi_alleff.py:52
tools.zlumi_alleff.trig_tag_and_probe
def trig_tag_and_probe(h, lb_minus_one_trig_eff)
Definition: zlumi_alleff.py:6
plotBeamSpotVxVal.range
range
Definition: plotBeamSpotVxVal.py:195
tools.zlumi_alleff.reco_tag_and_probe
def reco_tag_and_probe(hmo, hms, hno, hns, lb_minus_one_reco_eff, MC=False)
Definition: zlumi_alleff.py:19
dbg::print
void print(std::FILE *stream, std::format_string< Args... > fmt, Args &&... args)
Definition: SGImplSvc.cxx:70
tools.zlumi_alleff.extract
def extract(histogram, bin1, bin2)
Definition: zlumi_alleff.py:272
tools.zlumi_alleff.template_method_error
def template_method_error(hmo, hms, hno, hns, hto, hts)
Definition: zlumi_alleff.py:130
readCCLHist.float
float
Definition: readCCLHist.py:83
Generated on Thu Nov 7 2024 21:32:25 for ATLAS Offline Software by   doxygen 1.8.18