EvoMonitoring.py

Go to the documentation of this file.

# Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
#
#   This script is used to make plots comparing the latest IOV in the DB with the new calibration.
#   In order to run it standalone you need to setup Athena and call setupRunning(path_newCalib, path_oldCalib)
#   oldCalib file has the structure of the Recobery.py (after running the calibration)
#   newCalib file has the structure of the MakeReferenceFile (It is on IOV from the central DB)
#
 
from PixelCalibAlgs.Recovery import ReadDbFile, ReadNewCalib
from PathResolver import PathResolver
import array as ar
import numpy as np
from matplotlib.backends.backend_agg import FigureCanvasAgg
from matplotlib.figure import Figure
 
MaxThreshold = 10
MaxSlope = 10
 
def arrayCharge(parameters, layer):
    
    if layer == "Blayer":
        # Range [3, 6, 9, 12, 15, 18, 21, 24, 27, 30, 33, 36, 39, 42] (includes the ToT tuning point: 18)
        m_array = [charge(parameters, 3*(1+i)) for i in range(14)]
    else:
        # Range [5, 10, 15, 20, 25, 30, 35, 40] (includes the ToT tuning point: 30)
        m_array = [charge(parameters, 5*(1+i)) for i in range(8)]
        
    return ar.array('d',m_array)
 
def charge(parameters, tot):
    num = tot*parameters[2] - parameters[0]*parameters[1]
    den = parameters[0] - tot
    return (num/den) if den != 0 else 0
 
def percent( a,b):
    return a/(a+b)*100 if (a+b) != 0 else 999.99
 
 
def failsCheckTuning(parameters, layer):
    tot    = 30
    error  = MaxSlope # percentage
    expected_charge = 20000 # electrons
    
    if layer == "Blayer":
        tot = 18
        
    if layer == "IBL":
        tot = 10
        expected_charge = 16000 # electrons
        
        if abs(parameters[9]-expected_charge)/expected_charge*100 > 5:
            return True, [parameters[9], abs(parameters[9]-expected_charge)/expected_charge*100, expected_charge]
        return False, [0,0,0]
    
    low   = expected_charge*(1-error/100)
    high  = expected_charge*(1+error/100)
    realQ = charge(parameters, tot)
    
    if realQ < low or realQ > high:
        return True, [charge(parameters, tot), abs(charge(parameters, tot)-expected_charge)/expected_charge*100,expected_charge]
    return False, [0,0,0] 
 
def EvoMon(old_calib, new_calib, mapping, old_iov, new_iov):
    
    save = 1
    log_info = {}
    
    slopes = []
    information = {"Total_mods": 0, 
                   "Total_FE"  : 0, 
                   "IBL"   :{"bad":0, "ok":0}, 
                   "Blayer":{"bad":0, "ok":0}, 
                   "L1"    :{"bad":0, "ok":0}, 
                   "L2"    :{"bad":0, "ok":0}, 
                   "Disk"  :{"bad":0, "ok":0} 
                  }
    
    # There are 2048 modules in the structure
    for mod in range(2048):
 
        mod_str = mapping[str(mod)]
        mod_layer = ""
        print( "%-18s - %4i" % (mod_str, mod), end='\r')
        plot_range = np.array([0,75000])
        # Skipping modules that have not been calibrated
        if str(mod) not in new_calib:
            continue
               
        if mod_str.startswith("L0"): 
            mod_layer = "Blayer"
        elif mod_str.startswith("L1"): 
            mod_layer = "L1"
        elif mod_str.startswith("L2"): 
            mod_layer = "L2"
        elif mod_str.startswith("D"): 
            mod_layer = "Disk"
        else:
            mod_layer = "IBL"
            if mod_str.startswith("LI_S15"): 
                continue
        
        information["Total_mods"] += 1
        fig = Figure(figsize=(13,10))
        axs = fig.add_subplot(1,1,1)
        status = "_OK"
        
        for fe in range(len(new_calib[str(mod)])):
            
            information["Total_FE"] += 1
            newQ = []
            oldQ = []
            boolTOT = False
            realQ = []
            
            if mod_layer != "IBL":
                newCal_normal_pix = new_calib[str(mod)][fe][12:15]
                oldCal_normal_pix = old_calib[str(mod)][fe][12:15]
                boolTOT, realQ = failsCheckTuning(newCal_normal_pix, mod_layer)
            
                # We just fit the first point since we loose linearity afterwards
                newQ = arrayCharge(newCal_normal_pix,mod_layer)
                oldQ = arrayCharge(oldCal_normal_pix,mod_layer)
            else:
                # For IBL we dont need to convert TOT into charge, DB already in charge
                newQ = new_calib[str(mod)][fe][4:20]
                oldQ = old_calib[str(mod)][fe][4:20]
                
                plot_range = np.array([0,35000])
                boolTOT, realQ = failsCheckTuning(newQ, mod_layer)
                
            m,b = np.polyfit(newQ ,oldQ,1)
            slopes.append(m)
            boolFit = (abs((1-m)/m)*100) > MaxThreshold
            
            if boolFit or boolTOT:
                key = "%-18s - %i" % (mod_str, mod)
                if boolFit:
                    if key not in log_info:
                        log_info[key]  = "\tFE%02i ---> slope: %5.2f -  dev: %5.1f%%\n" % (fe,m, abs((1-m)/m)*100)
                    else:
                        log_info[key] += "\tFE%02i ---> slope: %5.2f -  dev: %5.1f%%\n" % (fe,m, abs((1-m)/m)*100)
                if boolTOT:
                    if key not in log_info:
                        log_info[key]  = "\tFE%02i ---> Charge= %5ie (dev %6.2f%%) out of error bars. Expected: %5ie\n" % (fe, realQ[0], realQ[1], realQ[2])
                    else:
                        log_info[key] += "\tFE%02i ---> Charge= %5ie (dev %6.2f%%) out of error bars. Expected: %5ie\n" % (fe, realQ[0], realQ[1], realQ[2])                    
                    
                information[mod_layer]["bad"] += 1
                status = "_BAD" 
                if boolTOT:
                    # Fails charge tuning
                    status += "_Q"
                if boolFit:
                    # Fails the slope at y = x (old vs. new calib)
                    status += "_Slope"
            else:
                information[mod_layer]["ok"] += 1
                
            
            if save:
                lstyle = "dotted" if fe < 8 else "solid"                                    
                axs.plot(newQ ,oldQ, marker='o', linestyle=lstyle, label = ("FE%02d" % (fe)))
                axs.set_xlim(plot_range)
                axs.set_ylim(plot_range)
            
        if save:
            fig.suptitle("Hash ID %d - %s" % (mod, mod_str))
            axs.legend(loc='upper left', ncol=4) 
            axs.set_xlabel(new_iov+" - Charge[e]")
            axs.set_ylabel(old_iov+" - Charge[e]")
            axs.plot(plot_range,plot_range,"k:",lw=1)
            
            if mod_layer == 'IBL':
                axs.plot([16000,16000],plot_range,"k:",lw=1)
                axs.text(17000,1000,"TOT@10 = 16 ke")
            elif mod_layer == 'Blayer':
                axs.plot([20000,20000],plot_range,"k:",lw=1)
                axs.text(21000,1000,"TOT@18 = 20 ke")
            else:
                axs.plot([20000,20000],plot_range,"k:",lw=1)
                axs.text(21000,1000,"TOT@30 = 20 ke")
            
            storage = "plots/" + mod_layer + "/"
            
            canvas = FigureCanvasAgg(fig)
            canvas.print_figure(storage+"Id"+str(mod)+ "_" +mod_str+status+".png", dpi=150)
    
    if(save):
        fig = Figure(figsize=(13,10))
        fig.suptitle("Fit slopes for all modules")
        axs = fig.add_subplot(1,1,1)
        axs.hist(np.clip(slopes, -1, 2), bins=100)
        axs.set_xlabel("Fit slope")
        axs.set_ylabel("Counts")
        # add text box for the statistics
        stats = (f'$\\mu$ = {np.mean(np.clip(slopes, -1, 2)):.3f}\n'
                 f'$\\sigma$ = {np.std(np.clip(slopes, -1, 2)):.3f}')
        bbox = dict(boxstyle='round', fc='blanchedalmond', ec='orange', alpha=0.5)
        axs.text(0.95, 0.07, stats, fontsize=9, bbox=bbox, transform=axs.transAxes)        
        
        FigureCanvasAgg(fig).print_figure("plots/FitSlopes.png", dpi=150)
 
    return information, log_info
 
    
def ReadCSV():
    mydict = {}
    with open(PathResolver.FindCalibFile("PixelCalibAlgs/mapping.csv")) as fp:
        lines = fp.readlines()
        for line in lines:
            arrayline = line.rstrip("\n").split(', ')
            mydict[arrayline[1]] = arrayline[0]
    return mydict
 
def ReadCalibOutput(file):
    mydict = {}
    with open(file) as fp:
        lines = fp.readlines()
        for line in lines:
            arrayline = line.rstrip("\n").split(' ')
            
            if arrayline[0] not in mydict:
                mydict[arrayline[0]] = []
            
            # Removing empty values. This is usually used for IBL calib file
            while("" in arrayline):
                arrayline.remove("")
            
            mydict[arrayline[0]].append([ int(arrayline[i]) if i<13 else float(arrayline[i]) for i in range(1,len(arrayline)) ])
 
    return mydict, "[latest,0]"
 
def setupRunEvo(path_newCalib, path_oldCalib):
    new_iov = "Latest IOV"
    old_iov = "Older IOV"
    
    print("Files chosen for the comparison:")
    
    print("New calibration: '%s'" % path_newCalib)
    new_calib, new_iov = "", ""
    if "PIX_FINAL_calibration_candidate.txt" in path_newCalib:
        new_calib, new_iov = ReadCalibOutput(path_newCalib)
    elif "ChargeCalib_" in path_newCalib:
        new_calib, new_iov = ReadCalibOutput(path_newCalib)
    else:
        new_calib, read_report = ReadNewCalib(path_newCalib)
    
    print("Old calibration: '%s'" % path_oldCalib)
    old_calib, old_iov = ReadDbFile(path_oldCalib)
 
    mapping = ReadCSV()
    
    import os
    os.makedirs("plots/Blayer", exist_ok=True)
    os.makedirs("plots/L1"    , exist_ok=True)
    os.makedirs("plots/L2"    , exist_ok=True)
    os.makedirs("plots/Disk"  , exist_ok=True)
    os.makedirs("plots/IBL"   , exist_ok=True)    
    
    information, log_info = EvoMon(old_calib, new_calib, mapping, old_iov, new_iov )    
 
    # writing log files
    fout = open("EvoMonitoring_log.txt", "w")
    
    fout.write("_______________ LOG _______________\n\n" )
    fout.write("Files chosen for the comparison:\n" )
    fout.write("New calibration: '%s'\n" % (path_newCalib) )
    fout.write("Old calibration: '%s'\n" % (path_oldCalib) )
    fout.write("Data base IOV:   %s\n\n" % (old_iov) )
    
    fout.write("-"*20+" SUMMARY "+"-"*20+"\n" )
    fout.write("%-20s: %5i\n"   % ("Total mods in det.", information["Total_mods"]))
    fout.write("%-20s: %5i\n"   % ("Total FE in det."  , information["Total_FE"]  ))
    fout.write("%-20s: %5i\n"   % ("Total FE IBL"      , information["IBL"]["ok"]   +information["IBL"]["bad"]   ))
    fout.write("%-20s: %5i\n"   % ("Total FE Blayer"   , information["Blayer"]["ok"]+information["Blayer"]["bad"]))
    fout.write("%-20s: %5i\n"   % ("Total FE L1"       , information["L1"]["ok"]    +information["L1"]["bad"]    ))
    fout.write("%-20s: %5i\n"   % ("Total FE L2"       , information["L2"]["ok"]    +information["L2"]["bad"]    ))
    fout.write("%-20s: %5i\n\n" % ("Total FE Disk"     , information["Disk"]["ok"]  +information["Disk"]["bad"]  ))
    
    fout.write(f'FrontEnds deviating more than {MaxThreshold}% of TOT vs charge gradient between new and previous calibration:\n')
    fout.write("%-11s: %-4i (%6.2f%%)\n"   % ("IBL FEs"   , information["IBL"]["bad"]   , percent(information["IBL"]["bad"]   ,information["IBL"]["ok"])   ))    
    fout.write("%-11s: %-4i (%6.2f%%)\n"   % ("Blayer FEs", information["Blayer"]["bad"], percent(information["Blayer"]["bad"],information["Blayer"]["ok"])))    
    fout.write("%-11s: %-4i (%6.2f%%)\n"   % ("L1 FEs"    , information["L1"]["bad"]    , percent(information["L1"]["bad"]    ,information["L1"]["ok"])    ))    
    fout.write("%-11s: %-4i (%6.2f%%)\n"   % ("L2 FEs"    , information["L2"]["bad"]    , percent(information["L2"]["bad"]    ,information["L2"]["ok"])    ))    
    fout.write("%-11s: %-4i (%6.2f%%)\n\n" % ("Disk FEs"  , information["Disk"]["bad"]  , percent(information["Disk"]["bad"]  ,information["Disk"]["ok"])  ))    
    
    fout.write("+"*20+" List of bad FE "+"+"*20+"\n" )
    fout.write("Expected TOT vs. charge for the different layers:\n")
    fout.write("%-10s: TOT@%2i = %2ike\n"   % ("IBL"       , 10, 16))
    fout.write("%-10s: TOT@%2i = %2ike\n"   % ("Blayer"    , 18, 20))
    fout.write("%-10s: TOT@%2i = %2ike\n\n" % ("L1/L2/Disk", 30, 20))
    
    for key, val in log_info.items():
        fout.write(key+"\n")
        fout.write(val+"\n")
        
    fout.close()
        
    #open and read the file for terminal print:
    fin = open("EvoMonitoring_log.txt", "r")
    print(fin.read())
 
 
if __name__ == "__main__":
 
    import argparse
    parser = argparse.ArgumentParser(prog='python -m PixelCalibAlgs.EvoMonitoring',
                            description="""Compares two IOV and plots the results.\n\n
                            Example: python -m PixelCalibAlgs.EvoMonitoring --new "path/to/file" --old "path/to/file" """)
    
    parser.add_argument('--new', default="PIX_FINAL_calibration_candidate.txt", help="New calibration file (output format from the Recovery.py)")
    parser.add_argument('--old', default="PixelChargeCalibration-DATA-RUN2-UPD4-28.log", help="Old DB IOV calibration")
    
    args = parser.parse_args()
    setupRunEvo(args.new, args.old)
    exit(0)