diff --git a/macros/leakageCurrent/getCurrent.py b/macros/leakageCurrent/getCurrent.py index cb33a28..9160c16 100755 --- a/macros/leakageCurrent/getCurrent.py +++ b/macros/leakageCurrent/getCurrent.py @@ -13,185 +13,185 @@ # obtaining for a given sector the integrated luminosity, the time and the measured current # as an array def vectors(sec): - time,lumi,curr = [],[],[] - for d in sec: - if d["current"]>0.00: - # excluding jumps in the measured currents by more than a factor 10 - if len(curr)>0 and d["current"]>10.0*curr[len(curr)-1]: - continue - if d["fill"] in bad_boys: - continue - lumi.append(d["intLumi"]) - time.append(d["endTime"]) - curr.append(d["current"]) + time,lumi,curr = [],[],[] + for d in sec: + if d["current"]>0.00: + # excluding jumps in the measured currents by more than a factor 10 + if len(curr)>0 and d["current"]>10.0*curr[len(curr)-1]: + continue + if d["fill"] in bad_boys: + continue + lumi.append(d["intLumi"]) + time.append(d["endTime"]) + curr.append(d["current"]) - return [lumi,time,curr] + return [lumi,time,curr] # function to correct for the temperature -> normalisation from a temperature t1 to T = 8 degrees C (t2) # Eg: silicon energy gap in eV def tempCorr(t1,t2=8.0,Eg=1.21): - kelvin = 273.15 - kb = 8.6173324e-5 - t1 += kelvin - t2 += kelvin - return (t1/t2)**2*np.exp(-Eg/(2*kb)*(1/t2-1/t1)) + kelvin = 273.15 + kb = 8.6173324e-5 + t1 += kelvin + t2 += kelvin + return (t1/t2)**2*np.exp(-Eg/(2*kb)*(1/t2-1/t1)) # function to apply the temperature correction for TT (i.e. two relevant temperature senssors CT5 and CT6 def corrCurrTT(time,curr,temp): - i1=0 - i2=0 - for i in range(0,len(curr)): - while time[i]>temp["CT5"][0][i1] and i1temp["CT6"][0][i2] and i2temp["CT5"][0][i1] and i1temp["CT6"][0][i2] and i2temp[0][j] and jtemp[0][j] and j=0) + return dict((key,value) for key, value in dic.iteritems() if key.find(expr)>=0) # filter high-voltage channels according to a given section/box name pattern def filterSect(dic,expr): - return dict((key,value) for key, value in dic.iteritems() if key.find(expr)==len(key)-len(expr)) + return dict((key,value) for key, value in dic.iteritems() if key.find(expr)==len(key)-len(expr)) # function to get the silicon volume connected to a particular value # returns 100000000000000.0 if the pattern does not match def getVolume(key): # in cm^3 - if key.find("TT")>=0: - oneSensor = 0.0500*9.64*9.44 - if key.find("US")>=0 or key.find("X1S")>=0: - if np.size(np.where(np.array([key.find(sec) for sec in ("S10","S4","S3","S9","S15","S16")])>=0))>0: - return 1.0*oneSensor - elif np.size(np.where(np.array([key.find(sec) for sec in ("S11","S5","S8","S14","S17")])>=0))>0: - return 2.0*oneSensor - elif np.size(np.where(np.array([key.find(sec) for sec in ("S12","S6","S7","S13","S18")])>=0))>0: - return 4.0*oneSensor - elif key.find("S2")==len(key)-2: - return 2.0*oneSensor - elif key.find("S1")==len(key)-2: - return 4.0*oneSensor - elif key.find("VS")>=0 or key.find("X2S")>=0: - if np.size(np.where(np.array([key.find(sec) for sec in ("S14","S8","S7","S13","S19","S20")])>=0))>0: - return 1.0*oneSensor - elif np.size(np.where(np.array([key.find(sec) for sec in ("S15","S9","S6","S12","S18","S21")])>=0))>0: - return 2.0*oneSensor - elif np.size(np.where(np.array([key.find(sec) for sec in ("S25","S24","S3")])>=0))>0: - return 3.0*oneSensor - elif np.size(np.where(np.array([key.find(sec) for sec in ("S26","S22","S16","S10","S4","S5","S11","S17","S23")])>=0))>0: - return 4.0*oneSensor - elif key.find("S2")==len(key)-2: - return 3.0*oneSensor - elif key.find("S1")==len(key)-2: - return 4.0*oneSensor + if key.find("TT")>=0: + oneSensor = 0.0500*9.64*9.44 + if key.find("US")>=0 or key.find("X1S")>=0: + if np.size(np.where(np.array([key.find(sec) for sec in ("S10","S4","S3","S9","S15","S16")])>=0))>0: + return 1.0*oneSensor + elif np.size(np.where(np.array([key.find(sec) for sec in ("S11","S5","S8","S14","S17")])>=0))>0: + return 2.0*oneSensor + elif np.size(np.where(np.array([key.find(sec) for sec in ("S12","S6","S7","S13","S18")])>=0))>0: + return 4.0*oneSensor + elif key.find("S2")==len(key)-2: + return 2.0*oneSensor + elif key.find("S1")==len(key)-2: + return 4.0*oneSensor + elif key.find("VS")>=0 or key.find("X2S")>=0: + if np.size(np.where(np.array([key.find(sec) for sec in ("S14","S8","S7","S13","S19","S20")])>=0))>0: + return 1.0*oneSensor + elif np.size(np.where(np.array([key.find(sec) for sec in ("S15","S9","S6","S12","S18","S21")])>=0))>0: + return 2.0*oneSensor + elif np.size(np.where(np.array([key.find(sec) for sec in ("S25","S24","S3")])>=0))>0: + return 3.0*oneSensor + elif np.size(np.where(np.array([key.find(sec) for sec in ("S26","S22","S16","S10","S4","S5","S11","S17","S23")])>=0))>0: + return 4.0*oneSensor + elif key.find("S2")==len(key)-2: + return 3.0*oneSensor + elif key.find("S1")==len(key)-2: + return 4.0*oneSensor + else: + if np.size(np.where(np.array([key.find(sec) for sec in ("P8","P4","P5","P1")])>=0))>0: + return 12.0*oneSensor + else: + return 9.0*oneSensor else: - if np.size(np.where(np.array([key.find(sec) for sec in ("P8","P4","P5","P1")])>=0))>0: - return 12.0*oneSensor - else: - return 9.0*oneSensor - else: - if key.find("BOXA")>=0 or key.find("BOXC")>=0: - return 0.041*7.6*11*2.0*4 - else: - return 0.032*7.6*11*4 - return 100000000000000.0 + if key.find("BOXA")>=0 or key.find("BOXC")>=0: + return 0.041*7.6*11*2.0*4 + else: + return 0.032*7.6*11*4 + return 100000000000000.0 # returns one multigraph for time and one multigraph for time for TT # filSector and filLayer are iterables representing different patterns to be selected # temp is a dictionnary of lists containing two elements, a time list and the corresponding temperature list def drawTT_sensor(data,filSector,filLayer,name,color,temp): - mg_time = ROOT.TMultiGraph("mg_%s_time"%name,"Multigraph time"); - mg_lumi = ROOT.TMultiGraph("mg_%s_lumi"%name,"Multigraph lumi"); - for f1 in filSector: - filter_dict = filterSect(data,f1) - for f2 in filLayer: - dfilter_dict = filterLayer(filter_dict,f2) - for sec in dfilter_dict: - print sec,": ",getVolume(sec) - [lumi,time,curr] = vectors(filter_dict[sec]) - corrCurrTT(time,curr,temp) - v_curr = np.array(curr)/getVolume(sec) - g_time = ROOT.TGraph(len(time),np.array(time)*1.0,v_curr) - g_lumi = ROOT.TGraph(len(time),np.array(lumi)*1.0,v_curr) - g_time.SetLineColor(color); - g_lumi.SetLineColor(color); - mg_time.Add(g_time,"l") - mg_lumi.Add(g_lumi,"l") - return [mg_time,mg_lumi] + mg_time = ROOT.TMultiGraph("mg_%s_time"%name,"Multigraph time"); + mg_lumi = ROOT.TMultiGraph("mg_%s_lumi"%name,"Multigraph lumi"); + for f1 in filSector: + filter_dict = filterSect(data,f1) + for f2 in filLayer: + dfilter_dict = filterLayer(filter_dict,f2) + for sec in dfilter_dict: + print sec,": ",getVolume(sec) + [lumi,time,curr] = vectors(filter_dict[sec]) + corrCurrTT(time,curr,temp) + v_curr = np.array(curr)/getVolume(sec) + g_time = ROOT.TGraph(len(time),np.array(time)*1.0,v_curr) + g_lumi = ROOT.TGraph(len(time),np.array(lumi)*1.0,v_curr) + g_time.SetLineColor(color); + g_lumi.SetLineColor(color); + mg_time.Add(g_time,"l") + mg_lumi.Add(g_lumi,"l") + return [mg_time,mg_lumi] # returns one multigraph for time and one multigraph for time for IT # filSector and filLayer are iterables representing different patterns to be selected # temp is a dictionnary of lists containing two elements, a time list and the corresponding temperature list def drawIT_sensor(data,filBox,name,color,temp): - mg_time = ROOT.TMultiGraph("mg_%s_time"%name,"Multigraph time"); - mg_lumi = ROOT.TMultiGraph("mg_%s_lumi"%name,"Multigraph lumi"); - for f1 in filBox: - filter_dict = filterLayer(data,f1) - for sec in filter_dict: - print sec,": ",getVolume(sec) - [lumi,time,curr] = vectors(filter_dict[sec]) - corrCurrIT(time,curr,temp[sec[0:-3]]) - v_curr = np.array(curr)/getVolume(sec) - g_time = ROOT.TGraph(len(time),np.array(time)*1.0,v_curr) - g_lumi = ROOT.TGraph(len(time),np.array(lumi)*1.0,v_curr) - g_time.SetLineColor(color); - g_lumi.SetLineColor(color); - mg_time.Add(g_time,"l") - mg_lumi.Add(g_lumi,"l") - return [mg_time,mg_lumi] + mg_time = ROOT.TMultiGraph("mg_%s_time"%name,"Multigraph time"); + mg_lumi = ROOT.TMultiGraph("mg_%s_lumi"%name,"Multigraph lumi"); + for f1 in filBox: + filter_dict = filterLayer(data,f1) + for sec in filter_dict: + print sec,": ",getVolume(sec) + [lumi,time,curr] = vectors(filter_dict[sec]) + corrCurrIT(time,curr,temp[sec[0:-3]]) + v_curr = np.array(curr)/getVolume(sec) + g_time = ROOT.TGraph(len(time),np.array(time)*1.0,v_curr) + g_lumi = ROOT.TGraph(len(time),np.array(lumi)*1.0,v_curr) + g_time.SetLineColor(color); + g_lumi.SetLineColor(color); + mg_time.Add(g_time,"l") + mg_lumi.Add(g_lumi,"l") + return [mg_time,mg_lumi] # combines an iterable of graphs into a multigraph def combineMultigraph(name,graphs): - mg = ROOT.TMultiGraph("mg_%s"%name,"Multigraph"); - for g in graphs: - mg.Add(g) - return mg + mg = ROOT.TMultiGraph("mg_%s"%name,"Multigraph"); + for g in graphs: + mg.Add(g) + return mg # sets the axis ranges etc. for a graph showing the current as a function of time def modifyGraphTime(mg): - #mg.GetXaxis().SetLimits(1262304000,1360000000); - mg.GetXaxis().SetLimits(1262304000,1438003626); - mg.GetXaxis().SetTimeDisplay(1); - mg.GetXaxis().SetTimeOffset(0.0); - mg.GetXaxis().SetTimeFormat("#splitline{%d/%m}{%Y}"); - mg.GetXaxis().SetTitle("Date"); - mg.GetYaxis().SetTitle("#it{I}_{leak}/#it{V} [A/m^{3}]"); - mg.GetXaxis().SetLabelOffset(0.02); - mg.GetXaxis().SetTimeOffset(0.0); - mg.GetYaxis().SetTitleOffset(1.3); - mg.GetXaxis().SetTitleOffset(1.4); + #mg.GetXaxis().SetLimits(1262304000,1360000000); + mg.GetXaxis().SetLimits(1262304000,1438003626); + mg.GetXaxis().SetTimeDisplay(1); + mg.GetXaxis().SetTimeOffset(0.0); + mg.GetXaxis().SetTimeFormat("#splitline{%d/%m}{%Y}"); + mg.GetXaxis().SetTitle("Date"); + mg.GetYaxis().SetTitle("#it{I}_{leak}/#it{V} [A/m^{3}]"); + mg.GetXaxis().SetLabelOffset(0.02); + mg.GetXaxis().SetTimeOffset(0.0); + mg.GetYaxis().SetTitleOffset(1.3); + mg.GetXaxis().SetTitleOffset(1.4); # sets the axis ranges etc. for a graph showing the current as a function of the luminosity def modifyGraphLumi(mg): - mg.GetXaxis().SetLimits(0.0,3500.0); - mg.GetXaxis().SetTitle("Delivered Luminosity [pb^{-1}]"); - mg.GetYaxis().SetTitle("#it{I}_{leak}/#it{V} [A/m^{3}]"); - mg.GetXaxis().SetLabelOffset(0.02); - mg.GetXaxis().SetTimeOffset(0.0); - mg.GetYaxis().SetTitleOffset(1.3); - mg.GetXaxis().SetTitleOffset(1.2); + mg.GetXaxis().SetLimits(0.0,3500.0); + mg.GetXaxis().SetTitle("Delivered Luminosity [pb^{-1}]"); + mg.GetYaxis().SetTitle("#it{I}_{leak}/#it{V} [A/m^{3}]"); + mg.GetXaxis().SetLabelOffset(0.02); + mg.GetXaxis().SetTimeOffset(0.0); + mg.GetYaxis().SetTitleOffset(1.3); + mg.GetXaxis().SetTitleOffset(1.2); # Running the lhcb style file