ATLAS Offline Software
LumiCalibrator.py
Go to the documentation of this file.
1 # Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
2 
3 #
4 # LumiCalibrator
5 #
6 # Eric Torrence - September 2011
7 #
8 # Contents:
9 # LumiCalibrator - utility class to apply the luminosity calibration stored in COOL
10 # to raw per-bcid lumi values.
11 #
12 # This must be initialized with the payload of a calibration record from COOL, then the raw luminosity
13 # will be calibrated with a call to calibrate
14 
15 import math
16 
17 from CoolLumiUtilities.LumiBlobConversion import bConvertList
18 
20 
21  def __init__(self):
22 
23  self.nPar = 0
24  self.fType = 'None'
25  self.muToLumi = 0.
26  self.parVec = []
27  self.rflag = 0
28 
29  self.currVersion = 0
30  self.currParA = 0.
31  self.currParB = 0.
32  self.currParC = 0.
33  self.currType = 1 # Use A+C by default
34 
35  self.polyVersion = 0
36  self.polyParA = 0.
37  self.polyParB = 0.
38  self.polyParC = 0.
39 
40  self.nBunches = 1
41 
42  # PMT currents
43  self.currentA = -1.
44  self.currentC = -1.
45 
46  self.oldCurrentA = -1.
47  self.oldCurrentC = -1.
48 
49  self.verbose = False
50  # self.verbose = True
51 
52  self.facTable = [1, 1]
53 
54  self.nParUsed = 0
55 
56  # Cached factorial values
57  def factorial(self, n):
58 
59  if n >= len(self.facTable):
60 
61  last = len(self.facTable) - 1
62  total = self.facTable[last]
63  for i in range(last+1, n+1):
64  total *= i
65  self.facTable.append(total)
66 
67  return self.facTable[n]
68 
69  def setCalibration(self, calibPayload):
70  self.nPar = calibPayload['NumOfParameters']
71  self.fType = calibPayload['Function']
72  self.muToLumi = calibPayload['MuToLumi']
73  blob = calibPayload['Parameters']
74  s=blob.read()
75  self.parVec = bConvertList(s, 4, self.nPar)
76 
77  if self.verbose:
78  print('LumiCalibrator.setCalibration for %s found %d parameters: ' % (self.fType, self.nPar))
79  print(self.parVec)
80  print(' muToLumi = %f' % self.muToLumi)
81 
82 
83  algname = self.fType
84  polyPos = algname.find("_PolyCorr")
85  currPos = algname.find("_CurrCorr")
86 
87  if polyPos > 0 and currPos > 0:
88  algname = algname[:-18]
89  elif polyPos > 0 or currPos > 0:
90  algname = algname[:-9]
91 
92  if self.verbose:
93  print('LumiCalibrator.setCalibration - Found poly: %d curr: %d' % (polyPos, currPos))
94 
95  if polyPos > currPos:
96 
97  if self.verbose:
98  print('LumiCalibrator.setCalibration - reading polynomial from', self.parVec[-4:-1])
99 
100  self.polyParC = self.parVec.pop()
101  self.polyParB = self.parVec.pop()
102  self.polyParA = self.parVec.pop()
103  self.polyVersion = self.parVec.pop()
104  self.nPar -= 4
105 
106  if currPos > 0:
107 
108  # Call calibrate to get number of parameters
109  self.calibrate(0.1)
110  npar = self.nParUsed
111  self.currType = 1 # Use A+C by default
112 
113  self.currVersion = int(self.parVec[npar])
114  if self.verbose:
115  print('LumiCalibrator.setCalibration - found %d parameters in %s' % (npar, algname))
116  print('LumiCalibrator.setCalibration - reading current vers %d from position %d' % (self.currVersion, npar))
117 
118  if self.currVersion == 4:
119  self.currType = int(self.parVec.pop())
120  self.nPar -= 1
121 
122  if self.currVersion >= 3:
123  self.currParC = self.parVec.pop()
124  self.nPar -= 1
125 
126  self.currParB = self.parVec.pop()
127  self.currParA = self.parVec.pop()
128  self.currVersion = self.parVec.pop()
129  self.nPar -= 3
130 
131  if self.verbose:
132  print('LumiCalibrator.setCalibration found current correction %d type %d' % (self.currVersion, self.currType))
133 
134  elif currPos > polyPos:
135 
136  # Call calibrate to get number of parameters
137  self.calibrate(0.1)
138  npar = self.nParUsed
139  self.currType = 1 # Use A+C by default
140 
141  if self.verbose:
142  print('LumiCalibrator.setCalibration - found %d parameters in %s' % (npar, algname))
143  if polyPos > 0:
144  npar += 4
145 
146  self.currVersion = int(self.parVec[npar])
147 
148  if self.verbose:
149  print('LumiCalibrator.setCalibration - reading current vers %d from position %d' % (self.currVersion, npar))
150 
151  if self.currVersion == 4:
152  self.currType = int(self.parVec.pop())
153  self.nPar -= 1
154 
155  if self.currVersion >= 3:
156  self.currParC = self.parVec.pop()
157  self.nPar -= 1
158 
159  self.currParB = self.parVec.pop()
160  self.currParA = self.parVec.pop()
161  self.currVersion = self.parVec.pop()
162  self.nPar -= 3
163 
164  if self.verbose:
165  print('LumiCalibrator.setCalibration found current correction %d type %d' % (self.currVersion, self.currType))
166 
167  if polyPos > 0:
168  if self.verbose:
169  print('LumiCalibrator.setCalibration - reading polynomial from list' , self.parVec[-4:])
170 
171  self.polyParC = self.parVec.pop()
172  self.polyParB = self.parVec.pop()
173  self.polyParA = self.parVec.pop()
174  self.polyVersion = self.parVec.pop()
175  self.nPar -= 4
176 
177 
178  def setLucidCurrent(self, currPayload):
179 
180  self.oldCurrentA = self.currentA
181  self.oldCurrentC = self.currentC
182 
183  if currPayload is None:
184  self.currentA = -1.
185  self.currentC = -1.
186  else:
187  self.currentA = currPayload['CurrentSideA']
188  self.currentC = currPayload['CurrentSideC']
189 
190 
191  if self.verbose:
192  print('LumiCalibrator.setLucidCurrent found Side A: %f Side C: %f' % (self.currentA, self.currentC))
193 
194  def dump(self):
195 
196  print('Calibration: %s, muToLumi: %f', (self.fType, self.muToLumi))
197  if self.fType == 'Logarithm':
198  print('par[0]: %f', self.parVec[0])
199 
200  elif self.fType == 'LookupTable_EventAND_Lin':
201  print('par[0]: %f', self.parVec[0])
202  print('par[1]: %f', self.parVec[1])
203  print('par[5]: %f', self.parVec[5])
204 
205  # Take raw luminosity (rate) value and return calibrated luminosity
206  def calibrate(self, rawLumi):
207 
208  cal = -1.
209  self.rflag = 0 # Success unless we fail
210 
211  if (rawLumi < 0.):
212  print('LumiCalibrator.calibrate(%f) - non-physical value!' % rawLumi)
213  self.rflag = 1
214  return 0.
215 
216  elif (rawLumi == 0.):
217  return 0.
218 
219  # Check for polynomial correction
220  calibstr = self.fType
221 
222  if calibstr.find('Polynomial') == 0:
223 
224  cal = self.calibPolynomial(rawLumi)
225 
226  elif calibstr.find('Logarithm') == 0:
227 
228  cal = self.calibLogarithm(rawLumi)
229 
230  elif calibstr.find('HitLogarithm') == 0:
231 
232  cal = self.calibHitLogarithm(rawLumi)
233 
234  elif calibstr.find('LookupTable_EventAND_Lin') == 0:
235 
236  try:
237  cal = self.calibLookupTable(rawLumi)
238  except Exception as e:
239  cal = 0.
240  self.rflag = 1
241  print('LumiCalibrator.calibLookupTable(%f) - Error: %s' % (rawLumi, e))
242 
243  elif calibstr.find('LookupTable_EventAND_Log') == 0:
244 
245  try:
246  cal = self.calibLookupTableLog(rawLumi)
247  except Exception as e:
248  cal = 0.
249  self.rflag = 1
250  print('LumiCalibrator.calibLookupTableLog(%f) - Error: %s' % (rawLumi, e))
251 
252  elif calibstr.find('LookupTable_EventANDFull_Log') == 0:
253 
254  try:
255  cal = self.calibLookupTableFullLog(rawLumi)
256  except Exception as e:
257  cal = 0.
258  self.rflag = 1
259  print('LumiCalibrator.calibLookupTableFullLog(%f) - Error: %s' % (rawLumi, e))
260 
261  elif calibstr.find('LookupTablePoisson_Lin') == 0:
262 
263  try:
264  cal = self.calibLookupTablePoisson(rawLumi)
265  except Exception as e:
266  cal = 0.
267  self.rflag = 1
268  print('LumiCalibrator.calibLookupTablePoisson(%f) - Error: %s' % (rawLumi, e))
269 
270  elif calibstr.find('LookupTableZeroPoisson_Lin') == 0:
271 
272  try:
273  cal = self.calibLookupTableZeroPoisson(rawLumi)
274  except Exception as e:
275  cal = 0.
276  self.rflag = 1
277  print('LumiCalibrator.calibLookupTableZeroPoisson(%f) - Error: %s' % (rawLumi, e))
278 
279  else:
280  print('LumiCalibrator.calibrate(%f) - Unknown calibration type %s!' % (rawLumi, self.fType))
281  self.rflag = 1
282  cal = 0.
283 
284  if self.rflag == 1:
285  return cal
286 
287  # Try to do polynomial and current correction in proper order
288  polyPos = calibstr.find("_PolyCorr")
289  currPos = calibstr.find("_CurrCorr")
290 
291  # No extra corrections
292  if polyPos == -1 and currPos == -1:
293  pass
294 
295  elif polyPos == -1 and currPos > 0:
296  cal = self.currentCorrection(cal)
297 
298  elif polyPos > 0 and currPos == -1:
299  cal = self.polyCorrection(cal)
300 
301  elif polyPos < currPos:
302  cal = self.polyCorrection(cal)
303  cal = self.currentCorrection(cal)
304 
305  elif currPos < polyPos:
306  cal = self.currentCorrection(cal)
307  cal = self.polyCorrection(cal)
308 
309  else:
310  print('LumiCalibrator.calibrate() - I am so confused: %s' % calibstr)
311 
312  return self.muToLumi * cal
313 
314  def currentCorrection(self, mu):
315 
316  correction = 1.
317 
318  current = -1.
319  if self.currType == 1:
320  current = self.currentA+self.currentC
321 
322  elif self.currType == 2:
323  current = self.currentA
324 
325  elif self.currType == 3:
326  current = self.currentC
327 
328  else:
329  print('LumiCalibrator.currentCorrection() - unknown current type: %d' % self.currType)
330  self.rflag = 1
331  return mu
332 
333  if self.currVersion == 0:
334  return mu
335 
336  if current == -1.:
337  if self.verbose:
338  print('LumiCalibrator.currentCorrection() - invalid LUCID currents found - A: %f C: %f' % (self.currentA, self.currentC))
339  print('LumiCalibrator.currentCorrection() - try using last values found - A: %f C: %f' % (self.oldCurrentA, self.oldCurrentC))
340 
341  self.rflag = 1 # Flag this as bad
342 
343  # Try using previous version
344  if self.currType == 1:
345  current = self.oldCurrentA+self.oldCurrentC
346 
347  elif self.currType == 2:
348  current = self.oldCurrentA
349 
350  elif self.currType == 3:
351  current = self.oldCurrentC
352 
353  if current == -1.:
354  return mu
355 
356  elif self.currVersion == 1:
357  correction = 100./(100.+self.currParA+self.currParB*current/self.nBunches)
358 
359  elif self.currVersion == 2:
360  correction = 100./(100.+self.currParA+self.currParB*current)
361 
362  elif self.currVersion == 3:
363  correction = 100./(100.+self.currParA+self.currParB*current+self.currParC*current**2)
364 
365  elif self.currVersion == 4: # Same as 3
366  correction = 100./(100.+self.currParA+self.currParB*current+self.currParC*current**2)
367 
368  else:
369  print('LumiCalibrator.currentCorrection() - unknown calibration version %d' % self.currVersion)
370 
371  if self.verbose and self.currVersion != 0:
372  print('LumiCalibrator.currentCorrection() - version %d -> currentA = %f, currentC = %f, correction = %f' % (self.currVersion, self.currentA, self.currentC, correction))
373 
374  return mu*correction
375 
376  def polyCorrection(self, mu):
377 
378  mucorr = mu
379 
380  # Parameters have already been parsed
381 
382  if self.polyVersion == 0:
383  pass
384 
385  elif self.polyVersion == 1:
386  mucorr = self.polyParA + self.polyParB * mu + self.polyParC * mu * mu
387 
388  else:
389  print('LumiCalibrator.polyCorrection() - unknown calibration version %d' % self.polyVersion)
390 
391  return mucorr
392 
393  def calibPolynomial(self, rawLumi):
394  cal = 0.
395  self.rflag = 1
396 
397  nrange = int(self.parVec[0])
398  self.nParUsed = 8 * nrange + 1
399 
400  for i in range(nrange):
401  rmin = self.parVec[i+1]
402  rmax = self.parVec[i+2]
403  if rawLumi < rmax and rawLumi >= rmin:
404  a = [self.parVec[i+3]]
405  a.append(self.parVec[i+4])
406  a.append(self.parVec[i+5])
407  a.append(self.parVec[i+6])
408  a.append(self.parVec[i+7])
409  a.append(self.parVec[i+8])
410  for k in range(6):
411  cal += a[k]*pow(rawLumi, k)
412 
413  self.rflag = 0
414  break
415 
416  if self.rflag == 1:
417  print('LumiCalibrator.calibPolynomial(%f) - Value out of range' % rawLumi)
418 
419  return cal
420 
421  def calibLogarithm(self, rawLumi):
422 
423  cal = 0.
424  self.rflag = 0
425 
426  self.nParUsed = 1
427 
428  invEff = self.parVec[0]
429 
430  try:
431  cal = -invEff*math.log(1.-rawLumi)
432 
433  except Exception:
434  cal = 0.
435  self.rflag = 1
436  # Don't print for simple saturation
437  if rawLumi != 1.:
438  print('LumiCalibrator.calibLogarithm(%f) - Unphysical input!' % rawLumi)
439 
440  return cal
441 
442  def calibHitLogarithm(self, rawLumi):
443 
444  cal = 0.
445 
446  invEff = self.parVec[0]
447  channels = self.parVec[1]
448  offset = self.parVec[2]
449  maxRawLumiperBX = self.parVec[3]
450 
451  self.nParUsed = 4
452 
453  # rawLumi can be > 1 here, check value against maxRawLumiPerBX
454  if rawLumi > maxRawLumiperBX:
455  self.rflag = 1
456  print('LumiCalibrator.calibHitLogarithm(%f) - input greater than max range %f!' % (rawLumi, maxRawLumiperBX))
457  return cal
458 
459  try:
460  cal = -invEff*math.log(1.-rawLumi/channels)/(1-offset)
461 
462  except Exception:
463  cal = 0.
464  self.rflag = 1
465  print('LumiCalibrator.calibHitLogarithm(%f) - Unphysical input!' % rawLumi)
466 
467  return cal
468 
469  # LookupTable_EventAND_Lin
470  def calibLookupTable(self, rawLumi):
471 
472  cal = 0.
473 
474  self.nParUsed = 6
475 
476  if rawLumi < 0.:
477  self.rflag = 1
478  print('LumiCalibrator.calibLookupTable(%f) - Unphysical input to LUT!' % rawLumi)
479 
480  elif rawLumi > 0.:
481  sigo = self.parVec[0]
482  siga = self.parVec[1]
483  # Try fast algorithm first
484  cal = self.parVec[5]*self.getMuvis(rawLumi, sigo, siga)
485  if self.rflag == 1:
486  # Try again with slower algorithm
487  cal = self.parVec[5]*self.getMuvis2(rawLumi, sigo, siga)
488 
489  return cal
490 
491  # LookupTable_EventAND_Log
492  def calibLookupTableLog(self, rawLumi):
493 
494  cal = 0.
495 
496  self.nParUsed = 8
497 
498  if (rawLumi < 0.) or (rawLumi >= 1.):
499  self.rflag = 1
500  print('LumiCalibrator.calibLookupTableLog(%f) - Unphysical input to LUT!' % rawLumi)
501 
502  elif rawLumi > 0.:
503  sigo = self.parVec[0]
504  siga = self.parVec[1]
505  # Try fast algorithm first
506  cal = self.parVec[7]*self.getMuvis(rawLumi, sigo, siga)
507  if self.rflag == 1:
508  # Try again with slower algorithm
509  cal = self.parVec[7]*self.getMuvis2(rawLumi, sigo, siga)
510 
511  return cal
512 
513  # From Vincent
514  def getMuvis(self, rawPerBX, sigo, siga):
515 
516  # Guess?
517  munew = 0.01
518 
519  b = sigo/siga
520  a = (b + 1)/2.
521 
522  # Set error, and clear below if we converge
523  self.iflag = 1
524 
525  # Set a fixed number of cycles, but break if we converge faster
526  for i in range(30):
527  mu = munew
528  y = 1 - 2 * math.exp(-a * mu) + math.exp(-b * mu)
529  dy = 2 * a * math.exp(-a * mu) - b * math.exp(-b * mu)
530 
531  munew = mu - (y-rawPerBX)/dy
532 
533  #print "Iteration: %d, Munew: %f, Muold: %f" % (i, munew, mu)
534  if munew <= 0.:
535  return 0.
536 
537  if math.fabs(munew-mu)/munew < 1.e-5:
538  self.iflag = 0
539  return munew
540 
541 
542  print('LumiCalibrator.calibLookupTable(', rawPerBX, ') - did not converge (Vincent method)!')
543  return munew
544 
545  # From Mika
546  def getMuvis2(self, rawPerBX, sigo, siga):
547  muvl=1e-10
548  muvu=100
549  muvm=10
550  sr=sigo/siga
551  rbxl=self.rpbx(sr,muvl)
552  rbxu=self.rpbx(sr,muvu)
553  rbxm=self.rpbx(sr,muvm)
554 
555  # Set error and clear below if we converge
556  self.rflag = 1
557 
558  # Check that starting value is in the valid range
559  if rawPerBX < rbxl or rawPerBX > rbxu:
560  print('LumiCalibrator.calibLookupTable(', rawPerBX, ') - raw lumi value outside of LUT range', rbxl, 'to', rbxu, '!')
561  return 0.
562 
563  # Restrict to a fixed number of iterations
564  for i in range(50):
565  if rbxl<rawPerBX and rbxm>rawPerBX:
566  rbxu=rbxm
567  muvu=muvm
568  muvm=0.5*(muvu+muvl)
569  else:
570  rbxl=rbxm
571  muvl=muvm
572  muvm=0.5*(muvu+muvl)
573 
574  rbxm = self.rpbx(sr, muvm)
575 
576  if (muvu-muvl)/muvl < 1e-5:
577  self.rflag = 0
578  return muvm
579 
580  # Did not converge
581  print('LumiCalibrator.calibLookupTable(', rawPerBX, ') - did not converge (Mika method)!')
582  return muvm
583 
584  def rpbx(self, sr, muvis):
585  return 1 - 2*math.exp(-(1+sr)*muvis/2) + math.exp(-sr*muvis)
586 
587  # LookupTable_EventANDFull_Log
588  def calibLookupTableFullLog(self, rawLumi):
589 
590  cal = 0.
591 
592  self.nParUsed = 9
593 
594  if (rawLumi < 0.) or (rawLumi >= 1.):
595  self.rflag = 1
596  print('LumiCalibrator.calibLookupTableFullLog(%f) - Unphysical input to LUT!' % rawLumi)
597 
598  elif rawLumi > 0.:
599  sigA = self.parVec[0]
600  sigC = self.parVec[1]
601  sigAND = self.parVec[2]
602 
603  # Try fast algorithm first
604  cal = self.parVec[8]*self.getMuvisFull(rawLumi, sigA, sigC, sigAND)
605  if self.rflag == 1:
606  # Try again with slower algorithm
607  cal = self.parVec[8]*self.getMuvisFull2(rawLumi, sigA, sigC, sigAND)
608 
609  return cal
610 
611  # Full version
612  def getMuvisFull(self, rawPerBX, sigA, sigC, sigAND):
613 
614  # Guess?
615  munew = 0.01
616 
617  ra = sigA/sigAND
618  rc = sigC/sigAND
619  b = ra+rc-1 # sigOR/sigAND
620 
621  # Set error, and clear below if we converge
622  self.iflag = 1
623 
624  # Set a fixed number of cycles, but break if we converge faster
625  for i in range(30):
626  mu = munew
627  y = 1 - math.exp(-ra * mu) - math.exp(-rc * mu) + math.exp(-b * mu)
628  dy = ra * math.exp(-ra * mu) + rc * math.exp(-rc * mu) - b * math.exp(-b * mu)
629 
630  munew = mu - (y-rawPerBX)/dy
631 
632  #print "Iteration: %d, Munew: %f, Muold: %f" % (i, munew, mu)
633  if munew <= 0.:
634  return 0.
635 
636  if math.fabs(munew-mu)/munew < 1.e-5:
637  self.iflag = 0
638  return munew
639 
640 
641  print('LumiCalibrator.calibLookupTable(', rawPerBX, ') - did not converge (Vincent method)!')
642  return munew
643 
644  # From Mika
645  def getMuvisFull2(self, rawPerBX, sigA, sigC, sigAND):
646  muvl=1e-10
647  muvu=100
648  muvm=10
649 
650  ra = sigA/sigAND
651  rc = sigC/sigAND
652  sr=ra+rc-1 # = sigOR/sigAND
653  rbxl=self.rpbxFull(ra, rc, sr,muvl)
654  rbxu=self.rpbxFull(ra, rc, sr,muvu)
655  rbxm=self.rpbxFull(ra, rc, sr,muvm)
656 
657  # Set error and clear below if we converge
658  self.rflag = 1
659 
660  # Check that starting value is in the valid range
661  if rawPerBX < rbxl or rawPerBX > rbxu:
662  print('LumiCalibrator.calibLookupTable(', rawPerBX, ') - raw lumi value outside of LUT range', rbxl, 'to', rbxu, '!')
663  return 0.
664 
665  # Restrict to a fixed number of iterations
666  for i in range(50):
667  if rbxl<rawPerBX and rbxm>rawPerBX:
668  rbxu=rbxm
669  muvu=muvm
670  muvm=0.5*(muvu+muvl)
671  else:
672  rbxl=rbxm
673  muvl=muvm
674  muvm=0.5*(muvu+muvl)
675 
676  rbxm = self.rpbxFull(ra, rc, sr, muvm)
677 
678  if (muvu-muvl)/muvl < 1e-5:
679  self.rflag = 0
680  return muvm
681 
682  # Did not converge
683  print('LumiCalibrator.calibLookupTable(', rawPerBX, ') - did not converge (Mika method)!')
684  return muvm
685 
686  def rpbxFull(self, ra, rc, sr, muvis):
687  return 1 - math.exp(-ra*muvis) - math.exp(-rc*muvis) + math.exp(-sr*muvis)
688 
689  # LookupTablePoisson_Lin':
690  def calibLookupTablePoisson(self, rawLumi):
691 
692  offset = self.parVec[0]
693  maxRawLumiperBX = self.parVec[1]
694  #mu_max = self.parVec[2]
695  #points = self.parVec[3]
696  nRefs = int(self.parVec[4])
697  ref = self.parVec[5:]
698 
699  self.nParUsed = 5+nRefs
700 
701  cal = 0.
702  self.rflag = 1
703 
704  if rawLumi > maxRawLumiperBX:
705  print('LumiCalibrator.calibLookupTablePoisson(%f) - input greater than max range %f!' % (rawLumi, maxRawLumiperBX))
706  return cal
707 
708  # Iteratively solve
709  munew = 0.001
710 
711  # Set a fixed number of cycles, but break faster if we converge
712  for i in range(30):
713  mu = munew
714  # No constant term in hit counting
715  y = 0
716  dy = 0
717  # Now add additional terms in sum, divide out factor of e^(-mu)
718  # Use running product to calculate mu^n/n!
719  fact = 1
720  for j in range(nRefs):
721  n = j+1
722  #term = ref[j] * math.pow(mu, n) / self.factorial(n)
723  #term = ref[j] * math.pow(mu, n) / math.factorial(n)
724  fact *= (mu/n)
725  term = ref[j] * fact
726  y += term
727  dy += term*(n/mu-1)
728 
729  # Here we are computing the probability of a hit, so we compare to the hit rate
730  munew = mu - (y-rawLumi/math.exp(-mu))/dy # Remember to put back in e^(-mu)
731  cal = munew
732 
733  if munew <= 0.:
734  cal = 0.
735  print('LumiCalibrator.calibLookupTablePoisson(%f) - failed to converge (went negative)!'% (rawLumi))
736  return cal
737 
738  if math.fabs(munew-mu)/munew < 1.e-5:
739  self.rflag = 0
740  return cal/(1-offset)
741 
742 
743  print('LumiCalibrator.calibLookupTablePoisson(%f) - failed to converge (trials)!'% (rawLumi))
744  return cal/(1-offset)
745 
746  # LookupTableZeroPoisson_Lin':
747  def calibLookupTableZeroPoisson(self, rawLumi):
748 
749  offset = self.parVec[0]
750  maxRawLumiperBX = self.parVec[1]
751  #mu_max = self.parVec[2]
752  #points = self.parVec[3]
753  nRefs = int(self.parVec[4])
754  ref = self.parVec[5:]
755 
756  self.nParUsed = 5+nRefs
757 
758  cal = 0.
759  self.rflag = 1
760 
761  if rawLumi > maxRawLumiperBX:
762  print('LumiCalibrator.calibLookupTableZeroPoisson(%f) - input greater than max range %f!' % (rawLumi, maxRawLumiperBX))
763  return cal
764 
765  # Iteratively solve
766  munew = 0.001
767 
768  # Set a fixed number of cycles, but break faster if we converge
769  for i in range(30):
770  mu = munew
771  # Constant term first (n=0) - divide out constant factor of e^(-mu)
772  y = 1
773  dy = -1
774  fact = 1
775  # Now add additional terms in sum
776  for j in range(nRefs):
777  n = j+1
778  fact *= (mu/n)
779  # term = ref[j] * math.pow(mu, n) / self.factorial(n)
780  term = ref[j] * fact
781  y += term
782  dy += term*(n/mu-1)
783 
784  # Here we are computing the probability of a *zero*, so compare to observed *zero* rate: (1-rawLumi)
785  zRate = (1-rawLumi)
786  munew = mu - (y-zRate/math.exp(-mu))/dy # Remember to put back in e^(-mu)
787  cal = munew
788 
789  if munew <= 0.:
790  print('LumiCalibrator.calibLookupTableZeroPoisson(%f) - failed to converge (negative)!'% (rawLumi))
791  cal = 0.
792  return cal
793 
794  if math.fabs(munew-mu)/munew < 1.e-5:
795  self.rflag = 0
796  return cal/(1-offset)
797 
798  print('LumiCalibrator.calibLookupTableZeroPoisson(%f) - failed to converge (trials)!'% (rawLumi))
799  return cal/(1-offset)
python.LumiCalibrator.LumiCalibrator.nBunches
nBunches
Definition: LumiCalibrator.py:40
python.LumiCalibrator.LumiCalibrator.dump
def dump(self)
Definition: LumiCalibrator.py:194
python.LumiCalibrator.LumiCalibrator.currentA
currentA
Definition: LumiCalibrator.py:43
python.LumiCalibrator.LumiCalibrator.calibrate
def calibrate(self, rawLumi)
Definition: LumiCalibrator.py:206
python.LumiCalibrator.LumiCalibrator.factorial
def factorial(self, n)
Definition: LumiCalibrator.py:57
python.LumiCalibrator.LumiCalibrator.currType
currType
Definition: LumiCalibrator.py:33
python.LumiCalibrator.LumiCalibrator.iflag
iflag
Definition: LumiCalibrator.py:523
python.LumiBlobConversion.bConvertList
def bConvertList(b, nbyte=1, nval=1)
Definition: LumiBlobConversion.py:37
python.LumiCalibrator.LumiCalibrator
Definition: LumiCalibrator.py:19
python.LumiCalibrator.LumiCalibrator.setCalibration
def setCalibration(self, calibPayload)
Definition: LumiCalibrator.py:69
python.LumiCalibrator.LumiCalibrator.currentCorrection
def currentCorrection(self, mu)
Definition: LumiCalibrator.py:314
python.LumiCalibrator.LumiCalibrator.calibLogarithm
def calibLogarithm(self, rawLumi)
Definition: LumiCalibrator.py:421
python.LumiCalibrator.LumiCalibrator.polyParA
polyParA
Definition: LumiCalibrator.py:36
python.LumiCalibrator.LumiCalibrator.polyParB
polyParB
Definition: LumiCalibrator.py:37
python.LumiCalibrator.LumiCalibrator.oldCurrentC
oldCurrentC
Definition: LumiCalibrator.py:47
python.LumiCalibrator.LumiCalibrator.calibLookupTablePoisson
def calibLookupTablePoisson(self, rawLumi)
Definition: LumiCalibrator.py:690
dumpHVPathFromNtuple.append
bool append
Definition: dumpHVPathFromNtuple.py:91
python.LumiCalibrator.LumiCalibrator.calibPolynomial
def calibPolynomial(self, rawLumi)
Definition: LumiCalibrator.py:393
python.LumiCalibrator.LumiCalibrator.setLucidCurrent
def setLucidCurrent(self, currPayload)
Definition: LumiCalibrator.py:178
python.LumiCalibrator.LumiCalibrator.__init__
def __init__(self)
Definition: LumiCalibrator.py:21
python.LumiCalibrator.LumiCalibrator.polyParC
polyParC
Definition: LumiCalibrator.py:38
python.LumiCalibrator.LumiCalibrator.oldCurrentA
oldCurrentA
Definition: LumiCalibrator.py:46
python.LumiCalibrator.LumiCalibrator.polyVersion
polyVersion
Definition: LumiCalibrator.py:35
python.LumiCalibrator.LumiCalibrator.parVec
parVec
Definition: LumiCalibrator.py:26
python.LumiCalibrator.LumiCalibrator.nParUsed
nParUsed
Definition: LumiCalibrator.py:54
python.LumiCalibrator.LumiCalibrator.verbose
verbose
Definition: LumiCalibrator.py:49
python.LumiCalibrator.LumiCalibrator.calibLookupTableFullLog
def calibLookupTableFullLog(self, rawLumi)
Definition: LumiCalibrator.py:588
python.LumiCalibrator.LumiCalibrator.polyCorrection
def polyCorrection(self, mu)
Definition: LumiCalibrator.py:376
python.LumiCalibrator.LumiCalibrator.rpbx
def rpbx(self, sr, muvis)
Definition: LumiCalibrator.py:584
python.LumiCalibrator.LumiCalibrator.currParC
currParC
Definition: LumiCalibrator.py:32
plotBeamSpotVxVal.range
range
Definition: plotBeamSpotVxVal.py:194
python.LumiCalibrator.LumiCalibrator.currentC
currentC
Definition: LumiCalibrator.py:44
python.LumiCalibrator.LumiCalibrator.calibLookupTableLog
def calibLookupTableLog(self, rawLumi)
Definition: LumiCalibrator.py:492
python.LumiCalibrator.LumiCalibrator.getMuvis
def getMuvis(self, rawPerBX, sigo, siga)
Definition: LumiCalibrator.py:514
python.LumiCalibrator.LumiCalibrator.getMuvisFull2
def getMuvisFull2(self, rawPerBX, sigA, sigC, sigAND)
Definition: LumiCalibrator.py:645
python.LumiCalibrator.LumiCalibrator.getMuvisFull
def getMuvisFull(self, rawPerBX, sigA, sigC, sigAND)
Definition: LumiCalibrator.py:612
print
void print(char *figname, TCanvas *c1)
Definition: TRTCalib_StrawStatusPlots.cxx:25
python.LumiCalibrator.LumiCalibrator.currParB
currParB
Definition: LumiCalibrator.py:31
python.LumiCalibrator.LumiCalibrator.currParA
currParA
Definition: LumiCalibrator.py:30
python.LumiCalibrator.LumiCalibrator.getMuvis2
def getMuvis2(self, rawPerBX, sigo, siga)
Definition: LumiCalibrator.py:546
python.LumiCalibrator.LumiCalibrator.muToLumi
muToLumi
Definition: LumiCalibrator.py:25
python.LumiCalibrator.LumiCalibrator.calibLookupTable
def calibLookupTable(self, rawLumi)
Definition: LumiCalibrator.py:470
python.LumiCalibrator.LumiCalibrator.currVersion
currVersion
Definition: LumiCalibrator.py:29
python.CaloAddPedShiftConfig.int
int
Definition: CaloAddPedShiftConfig.py:45
python.LumiCalibrator.LumiCalibrator.calibLookupTableZeroPoisson
def calibLookupTableZeroPoisson(self, rawLumi)
Definition: LumiCalibrator.py:747
python.LumiCalibrator.LumiCalibrator.rflag
rflag
Definition: LumiCalibrator.py:27
python.LumiCalibrator.LumiCalibrator.nPar
nPar
Definition: LumiCalibrator.py:23
python.LumiCalibrator.LumiCalibrator.fType
fType
Definition: LumiCalibrator.py:24
python.LumiCalibrator.LumiCalibrator.calibHitLogarithm
def calibHitLogarithm(self, rawLumi)
Definition: LumiCalibrator.py:442
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
python.LumiCalibrator.LumiCalibrator.rpbxFull
def rpbxFull(self, ra, rc, sr, muvis)
Definition: LumiCalibrator.py:686
python.LumiCalibrator.LumiCalibrator.facTable
facTable
Definition: LumiCalibrator.py:52