LArHVGainsPredictor.py

Go to the documentation of this file.

# Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
 
from __future__ import print_function
 
from optparse import OptionParser
 
import PlotCalibrationHV
 
import PlotCalibrationGains
 
from PyCool import cool
 
import time
 
import sys
 
import mergeEnergyRamps
 
import ROOT
 
import os
 
 
 
class OracleGainReader:
 
  def __init__(self):
 
    self.reference_gains={}
    self.UNIX2COOL = 1000000000
 
    # get database service and open database
    dbSvc = cool.DatabaseSvcFactory.databaseService()
 
    dbString = 'oracle://ATLAS_COOLPROD;schema=ATLAS_COOLONL_TRIGGER;dbname=CONDBR2'
    try:
      db = dbSvc.openDatabase(dbString, False)        
    except Exception as e:
      print ('Error: Problem opening database', e)
      sys.exit(1)
 
    folder_name = "/TRIGGER/Receivers/Factors/CalibGains"
    folder=db.getFolder(folder_name)
       
    startUtime = int(time.time())
    endUtime = int(time.time())
    startValKey = startUtime * self.UNIX2COOL
    endValKey = endUtime * self.UNIX2COOL
    chsel = cool.ChannelSelection(0,sys.maxint)
 
    try:
      itr=folder.browseObjects(startValKey, endValKey, chsel)
    except Exception as e:
      print (e)
      sys.exit(1)
 
    for row in itr:
      ReceiverId = hex(int(row.channelId()))
      payload = row.payload()
      gain = payload['factor']
      self.reference_gains[ReceiverId]=gain
 
    # close database
    db.closeDatabase()
 
 
  def getGain(self,ReceiverId):
           
    if (ReceiverId in self.reference_gains):
      return float(self.reference_gains[ReceiverId])
    else:
      return ''
 
 
class GainPredictor:
 
  def __init__(self, mapping):
    
    self.geometry_convertor = mapping
 
    self.layer_weights_em  = {} # per eta bin
    self.layer_weights_had = {} # per eta bin
 
    file_name= ROOT.PathResolver.find_calib_file("TrigT1Calo/HVcorrPhysicsWeights_v1.txt")
 
    
    
    try:
      
      file = open(file_name)
      
    except IOError:
 
      print ("\ncould not find file: %s ....exiting\n" % file_name)
 
      sys.exit()
 
    for line in file.readlines():
 
      parts = line.split()
 
      self.layer_weights_em [int(parts[0])] = [float(parts[1]),float(parts[2]),float(parts[3]),float(parts[4]),float(parts[5])]
      self.layer_weights_had[int(parts[0])] = [float(parts[6]),float(parts[7]),float(parts[8]),float(parts[9])]
 
    file.close()
    
 
  def GetGain(self, receiver, orig_gain, layer_corr_ref, layer_corr_new, layer_names):
 
    coolid = self.geometry_convertor.getPPMfromReceiver(receiver)
    
    eta_bin = self.geometry_convertor.getEtaBin(coolid)
 
    hv_coef_ref = 0. # ref hv coefficient
    hv_coef_new = 0. # new hv coefficient
 
    
 
    if self.geometry_convertor.isCoolEm (coolid):
 
      is_overlap = self.geometry_convertor.isPPMOverlap(coolid)
 
      
 
      for n in range(len(layer_names)):
 
        layer = layer_names[n]
        
        if PlotCalibrationHV.isEmPresampler(layer):
 
          hv_coef_ref += layer_corr_ref[n] * self.layer_weights_em [eta_bin][0]
          hv_coef_new += layer_corr_new[n] * self.layer_weights_em [eta_bin][0]
 
        if PlotCalibrationHV.isEmFront(layer,is_overlap):
 
          hv_coef_ref += layer_corr_ref[n] * self.layer_weights_em [eta_bin][1]
          hv_coef_new += layer_corr_new[n] * self.layer_weights_em [eta_bin][1]
 
        if PlotCalibrationHV.isEmMiddle(layer,is_overlap):
 
          hv_coef_ref += layer_corr_ref[n] * self.layer_weights_em [eta_bin][2]
          hv_coef_new += layer_corr_new[n] * self.layer_weights_em [eta_bin][2]
 
        if PlotCalibrationHV.isEmBack(layer,is_overlap):
 
          hv_coef_ref += layer_corr_ref[n] * self.layer_weights_em [eta_bin][3]
          hv_coef_new += layer_corr_new[n] * self.layer_weights_em [eta_bin][3]
 
        if PlotCalibrationHV.isEmOverlapBack(layer,is_overlap):
 
          hv_coef_ref += layer_corr_ref[n] * self.layer_weights_em [eta_bin][4]
          hv_coef_new += layer_corr_new[n] * self.layer_weights_em [eta_bin][4]          
 
    
   
    if self.geometry_convertor.isCoolHad(coolid):
 
      regions = "EmbFcalHighEta" # (?)
 
      
 
      for n in range(len(layer_names)):
 
        layer = layer_names[n]
 
        if PlotCalibrationHV.isHadFirstLayer(layer,regions):
        
          hv_coef_ref += layer_corr_ref[n] * self.layer_weights_had[eta_bin][0]
          hv_coef_new += layer_corr_new[n] * self.layer_weights_had[eta_bin][0]
 
        if PlotCalibrationHV.isHadSecondLayer(layer,regions):
            
          hv_coef_ref += layer_corr_ref[n] * self.layer_weights_had[eta_bin][1]
          hv_coef_new += layer_corr_new[n] * self.layer_weights_had[eta_bin][1]
          
        if PlotCalibrationHV.isHadThirdLayer(layer,regions):
 
          hv_coef_ref += layer_corr_ref[n] * self.layer_weights_had[eta_bin][2]
          hv_coef_new += layer_corr_new[n] * self.layer_weights_had[eta_bin][2]
          
        if PlotCalibrationHV.isHadFourthLayer(layer,regions):
 
          hv_coef_ref += layer_corr_ref[n] * self.layer_weights_had[eta_bin][3]
          hv_coef_new += layer_corr_new[n] * self.layer_weights_had[eta_bin][3]
 
    pred_gain = orig_gain * (hv_coef_new/hv_coef_ref) # predicted gain
 
    return pred_gain
 
 
if __name__ == "__main__":
   
  
 
  parser = OptionParser(add_help_option=False)
  
  parser.add_option("-n", "--hv_input_new", action = "store", type = "string", dest = "hv_input_new", default =             "") 
  parser.add_option("-r", "--hv_input_ref", action = "store", type = "string", dest = "hv_input_ref", default =             "") 
  parser.add_option("-t", "--hv_corr_diff", action = "store", type =  "float", dest = "hv_corr_diff", default =           0.01)
  parser.add_option("-c", "--channel_list", action = "store", type = "string", dest = "channel_list", default =             "")
  parser.add_option("-o", "--output_files", action = "store", type = "string", dest = "output_files", default = "gain_updates")
 
  parser.add_option("-h", "--help", action = "store_true")
 
  (options, args) = parser.parse_args()
 
  if options.help: # print helpful info 
    
    print ("\nusage:    python LArHVGainsPredictor.py [options]")
    
    print ("\noptions:\n"                  )
 
    print ("-n, --hv_input_new    new hv input .sqlite file (see below)    (default:           '') *")
    print ("-r, --hv_input_ref    ref hv input .sqlite file (see below)    (default:           '') *")
    print ("-t, --hv_corr_diff    minimum abs. change in hv corrections    (default:         0.01)")
    print ("-c, --channel_list    file containing an input channel list    (default:           '')")
    print ("-o, --output_files    name assigned to all the output files    (default: gains_update)")
 
    print ("\n[*] minimal requirements")
 
    print ("\nexample:  python LArHVGainsPredictor.py -n file_1 -r file_2")
 
    print ("\nnote: if no input channel list is given all the channels will be considered by default")
 
    print ("\n(see https://twiki.cern.ch/twiki/bin/save/Atlas/LevelOneCaloGainPredictor for more info)\n")
 
    sys.exit()
 
  error_code = 2 # which error code (?)
 
  
 
  if options.hv_input_new == "" \
  or options.hv_input_ref == "":
 
    print ("\ntoo few input arguments given ....exiting  (see python LArHVGainsPredictor.py -h for more info)\n")
        
    sys.exit()
    
  hv_input_new = PlotCalibrationHV.L1CaloHVReader(options.hv_input_new) 
 
  hv_input_ref = PlotCalibrationHV.L1CaloHVReader(options.hv_input_ref)
 
  print ("\nreading hv input: %s and hv reference: %s" % (options.hv_input_new, options.hv_input_ref))
  
  
  
  geometry_convertor = PlotCalibrationGains.L1CaloGeometryConvertor()
 
  geometry_convertor.LoadReceiverPPMMap()
 
  
  
  gain_predictor = GainPredictor(geometry_convertor)
 
  
 
  h_orig_gains_em  = PlotCalibrationGains.L1CaloMap("Em:  orig. gains"," #eta bin","#phi bin")
  h_orig_gains_had = PlotCalibrationGains.L1CaloMap("Had: orig. gains"," #eta bin","#phi bin")
  
  h_pred_gains_em  = PlotCalibrationGains.L1CaloMap("Em:  pred. gains"," #eta bin","#phi bin")
  h_pred_gains_had = PlotCalibrationGains.L1CaloMap("Had: pred. gains"," #eta bin","#phi bin")
  
  h_diff_gains_em  = PlotCalibrationGains.L1CaloMap("Em:  (pred. gain - orig. gain) / orig. gain","#eta bin","#phi bin")
  h_diff_gains_had = PlotCalibrationGains.L1CaloMap("Had: (pred. gain - orig. gain) / orig. gain","#eta bin","#phi bin")
 
  
  
  output_text = open(options.output_files+".txt","w") 
 
  
  
  pred_gains = {}
 
  
 
  if options.channel_list != "": # from .txt file
 
    try:
 
      channel_list = open(options.channel_list)
 
    except IOError:
 
      print ("\ncould not find file: %s ....exiting\n" % options.channel_list)
 
      sys.exit()
 
    print ("\nreading input channel list from:", options.channel_list)
   
    orig_gains = {}
 
    for channel in channel_list.readlines():
 
      parts = channel.split()
    
      orig_gains[parts[0]] = float(parts[1])
 
    channel_list.close()
    
    receiver_list = orig_gains.keys()
 
  else: # from oracle database
 
    print ("\nreading input channel list from oracle database")
    
    orig_gains = OracleGainReader()
    
    receiver_list = geometry_convertor.receiver_to_ppm_map.keys()
 
  
 
  if not receiver_list:
 
    print ("\ninput channel list is empty ....exiting\n")
 
    sys.exit()
    
  print ("\nsearching %s channels for hv correction changes > %s" % (len(receiver_list), options.hv_corr_diff))
 
  
    
  print ("\npreparing gain updates for the following list of channels:")
    
  for receiver in receiver_list:
    
    if receiver not in hv_input_ref.GetNLayers().keys():
 
      continue # skip this receiver (it's non-LAr)
 
    coolid = geometry_convertor.getPPMfromReceiver(receiver)
        
    eta_bin = geometry_convertor.getEtaBin(coolid)
    phi_bin = geometry_convertor.getPhiBin(coolid)
    
    
 
    num_layers = (hv_input_ref.GetNLayers())[receiver]
 
    layer_names = [-1,-1,-1,-1] # default layer names (see... Calorimeter/CaloIdentifier/CaloIdentifier/CaloCell_ID.h)
    
    if num_layers > 0:
      layer_names[0] = (hv_input_ref.GetName1()[receiver])
    if num_layers > 1:
      layer_names[1] = (hv_input_ref.GetName2()[receiver])
    if num_layers > 2:
      layer_names[2] = (hv_input_ref.GetName3()[receiver])
    if num_layers > 3:
      layer_names[3] = (hv_input_ref.GetName4()[receiver])
  
    
 
    layer_corr_new = [1.,1.,1.,1.] # default hv corrections
 
    if receiver in hv_input_new.GetMeanCorections().keys(): # if mean hv correction > 1
      
      if num_layers > 0:
        layer_corr_new[0] = (hv_input_new.GetCorLayer1()[receiver])
      if num_layers > 1:
        layer_corr_new[1] = (hv_input_new.GetCorLayer2()[receiver])
      if num_layers > 2:
        layer_corr_new[2] = (hv_input_new.GetCorLayer3()[receiver])
      if num_layers > 3:
        layer_corr_new[3] = (hv_input_new.GetCorLayer4()[receiver])
        
    
 
    layer_corr_ref = [1.,1.,1.,1.] # default hv corrections
 
    if receiver in hv_input_ref.GetMeanCorections().keys(): # if mean hv correction > 1
      
      if num_layers > 0:
        layer_corr_ref[0] = (hv_input_ref.GetCorLayer1()[receiver])
      if num_layers > 1:
        layer_corr_ref[1] = (hv_input_ref.GetCorLayer2()[receiver])
      if num_layers > 2:
        layer_corr_ref[2] = (hv_input_ref.GetCorLayer3()[receiver])
      if num_layers > 3:
        layer_corr_ref[3] = (hv_input_ref.GetCorLayer4()[receiver])
        
    
 
    if  abs(layer_corr_new[0] - layer_corr_ref[0]) <= options.hv_corr_diff \
    and abs(layer_corr_new[1] - layer_corr_ref[1]) <= options.hv_corr_diff \
    and abs(layer_corr_new[2] - layer_corr_ref[2]) <= options.hv_corr_diff \
    and abs(layer_corr_new[3] - layer_corr_ref[3]) <= options.hv_corr_diff:
 
      continue # skip this receiver
     
    
      
    if options.channel_list != "": # from .txt file
      
      orig_gain = orig_gains[receiver]
 
    else: # from oracle database
        
      orig_gain = orig_gains.getGain(receiver)
        
    
 
    pred_gain = gain_predictor.GetGain(receiver,orig_gain,layer_corr_ref,layer_corr_new,layer_names)
 
    if not pred_gains:
      
      print ("") # just some formatting
    
    pred_gains[receiver] = [pred_gain, error_code]
      
    diff_gain = ((pred_gain - orig_gain) / orig_gain) * 100.0 # % change in gain
 
    
      
    if geometry_convertor.isCoolEm (coolid):
 
      layer =  "EM"
  
      h_orig_gains_em. Fill(eta_bin,phi_bin,orig_gain)
      h_pred_gains_em. Fill(eta_bin,phi_bin,pred_gain)
      h_diff_gains_em. Fill(eta_bin,phi_bin,diff_gain)
          
    if geometry_convertor.isCoolHad(coolid):
      
      layer = "HAD" 
      
      h_orig_gains_had.Fill(eta_bin,phi_bin,orig_gain)
      h_pred_gains_had.Fill(eta_bin,phi_bin,pred_gain)
      h_diff_gains_had.Fill(eta_bin,phi_bin,diff_gain)
      
    
 
    print (("%9s   %3s   %3i %2i   %6s   %.3f   %.3f   %6.1f %%   %1i") % (coolid, layer, eta_bin, phi_bin, receiver, orig_gain, pred_gain, diff_gain, num_layers), end='')
 
    print (("  [%s]") % ", ".join("%.3f" % x for x in layer_corr_ref), end='')
    
    print (("  [%s]") % ", ".join("%.3f" % x for x in layer_corr_new), end='')
 
    print (("  [%s]") % ", ".join( "%2i" % x for x in layer_names))
 
    
    
    print (("%6s   %.3f   %s") % (receiver, pred_gain, " ".join("%.3f" % x for x in layer_corr_new)), file=output_text, end='')
 
    print (("  #  "), file=output_text, end='')
 
    print (("%9s   %3i %2i   %3s   %.3f") % (coolid, eta_bin, phi_bin, layer, orig_gain), file=output_text, end='')
 
    print (("  %s") % " ".join("%6.3f" % (layer_corr_new[x] - layer_corr_ref[x]) for x in range (4)), file=output_text)
 
  output_text.close()
 
  
 
  mergeEnergyRamps.WriteSqlite(options.output_files+".sqlite",pred_gains)
 
  
  
  print ("\nrecreating text and ps/pdf files: %s and %s" % (options.output_files+".txt", options.output_files+".ps/pdf"))
 
  canvas = ROOT.TCanvas("canvas","",200,10,700,500)
 
  ROOT.gStyle.SetPalette(1)
  ROOT.gStyle.SetOptStat(111111)
  ROOT.gStyle.SetCanvasColor(10)
 
  canvas.Print(options.output_files+".ps[")
  
  h_orig_gains_em.Draw()
  canvas.Print(options.output_files+".ps")
 
  h_orig_gains_had.Draw()
  canvas.Print(options.output_files+".ps")
  
  h_pred_gains_em.Draw()
  canvas.Print(options.output_files+".ps")  
   
  h_pred_gains_had.Draw()
  canvas.Print(options.output_files+".ps")  
 
  h_diff_gains_em.Draw()
  canvas.Print(options.output_files+".ps")
 
  h_diff_gains_had.Draw()
  canvas.Print(options.output_files+".ps")
 
  canvas.Print(options.output_files+".ps]")
  
  os.system("ps2pdf "+options.output_files+".ps")
  
  print ("\nfinished.... need a drink now!\n")