Skip to content
Snippets Groups Projects
Commit 31837c84 authored by Samuel Rodriguez Ochoa's avatar Samuel Rodriguez Ochoa
Browse files

Upload New File

parent c5c708fb
Branches
No related tags found
No related merge requests found
#!/usr/bin/env python
#imports
import ROOT
import argparse
import numpy as np
from datetime import datetime
#Run code
#python -i 2024.05.17_analysis_SCB1.5_total_board.py -d 2024.05.30_DAQ_report.txt -e Samuel_Rodriguez -b SCB1
#---------------------------------------------------FUNCTIONS------------------------------------------------------------------
#Temperature function
def Temperature(R_SCB):
A=3.9083E-3
B=-5.775E-7
C=-4.183E-12
R0=10000
if R_SCB < R0:
# 0 = (C4)T^4 + (C3)T^3 + (C2)T^2 + (C1)T + (C0)
C4 = C*R0
C3 = -C*R0*100
C2 = B*R0
C1 = A*R0
C0 = R0-R_SCB
Tempf = ROOT.Math.Polynomial(C4,C3,C2,C1,C0)
else:
# 0 = (C2)T^2 + (C1)T + (C0)
C2 = B*R0
C1 = A*R0
C0 = R0-R_SCB
Tempf = ROOT.Math.Polynomial(C2,C1,C0)
return Tempf
pass
#---------------------------------------------------VARIABLES-----------------------------------------------------------------
#MEASUREMENTS
Ch = 0 # channel read
Iout = 0 # source current from the VRB in uA
Vd = 0 # experimental Vref measured
Vout = 0 # output voltage from SCB
Vin = 0 # input voltage from main power supply
Vref=0 # theoretical Vref then is set
points=0 #points taken
#FLAG VARIABLES
current_Vin=0 # when Vin change
sample_counter=0 # samples registered
flag_data = False # flag in false
#error
maximum_error=0 #Maximum delta measured collector
current_e1=0#maximum error from VRB
current_e2=0#maximum error from SCB
channel_error1=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]#channel error map
channel_error2=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]#channel error map
#--------------------------------------------------STYLE IMPORT---------------------------------------------------------
#ATLAS style
ROOT.gROOT.SetStyle("ATLAS")
legend = ROOT.TLegend(.85,.60,.95,.95)
#Colors for every channel
color_map=[1,11,3,4,5,6,7,8,9,30,38,40,41,42,43,44,45,47,28,23,21,12,33,16]
MuxMap=[4,5,8,9,12,13,14,3,6,7,10,11,27,26,23,22,19,18,17,28,25,24,21,20]
#---------------------------------------------------GRAPHS----------------------------------------------------------------
#multigraph Graph format
mg1 = ROOT.TMultiGraph()# Voltage (SCB) VS Temperature (VRB)
mg2 = ROOT.TMultiGraph()# Voltage (SCB)VS Temperature (SCB)
mg3 = ROOT.TMultiGraph()# Temperature Error (VRB) VS Temperature (VRB)
mg4 = ROOT.TMultiGraph()# Temperature Error (SCB) VS Temperature (SCB)
mg5 = ROOT.TMultiGraph()# Delta Temperature (VRB-SCB) VS Temperature (VRB)
#multigraph names
mg1.SetTitle(";Voltage from SCB [V];Temperature Read from VRB [C]")# Temperature (VRB) VS Voltage (SCB)
mg2.SetTitle(";Voltage from SCB [V];Temperature Read from SCB [C]")# Temperature (SCB) VS Voltage (SCB)
mg3.SetTitle(";Theoretical temperature [C];Temperature Error Read from VRB [C]")# Temperature Error (VRB) VS Temperature (THEORETICAL)
mg4.SetTitle(";Theoretical temperature [C];Temperature Error Read from SCB [C]")# Temperature Error (SCB) VS Temperature (THEORETICAL)
mg5.SetTitle(";Input from VRB temperature [C];#Delta Temperature [C]")# Delta Temperature (VRB-SCB) VS Temperature (VRB)
#data management
data1={}#temp (VRB)
data2={}#temp (SCB)
data3={}#error temp (VRB)
data4={}#error temp (SCB)
data7={}#Delta temperature (VRB-SCB)
#tf to fit
tf1={}# Voltage (SCB) VS Temperature (VRB)
tf2={}# Voltage (SCB)VS Temperature (SCB)
#-------------------------------------------------HISTOGRAMS-----------------------------------------------------------------
#analysis Voltage (SCB) VS Temperature (VRB) and Voltage (SCB)VS Temperature (SCB) parameters
h1=ROOT.TH1F("h1",";P0_VRB [C]",100,1000,2000)
h2=ROOT.TH1F("h2",";P1_VRB [C/V]",100,-2000,-1000)
h3=ROOT.TH1F("h3",";P2_VRB [C/V^2]",100,300,700)
h4=ROOT.TH1F("h4",";P0_SCB [C]",100,1000,2000)
h5=ROOT.TH1F("h5",";P1_SCB [C/V]",100,-2000,-1000)
h6=ROOT.TH1F("h6",";P2_SCB [C/V^2]",100,300,700)
#Error Delta Temperature (VRB-SCB)
h9=ROOT.TH1F("Delta_Temperature",";#Delta Temperature [C]",50,-1,1)
#Vref
h10=ROOT.TH1F("Vref_histogram",";Vref [V]",10,0.7,0.9) #Vref
#------------------------------------------------------MAIN----------------------------------------------------------
#-------------------------------data file---------------------------
parser = argparse.ArgumentParser()
parser.add_argument("-d","--input",required=True)#data file read
parser.add_argument("-e","--examinator",required=True)#examinator
parser.add_argument("-b","--SCB",required=True)#examinator
args = parser.parse_args()
#--------data collector
file=open(args.input,"r")
for line in file.readlines():
data = line.strip().split(',')#split by commas line read
ch = int(data[0])#channel
Vin = float(data[1])#Vin (negative)
Iout += float(data[2])#Iout
Vout += float(data[3])#Vout
Vd += float(data[4])#Vd
Vref= 0.799319 #Vref theoretical
sample_counter+=1
#--------data processor
#-------mean in every measure (voltage point)
if flag_data:
Iout=Iout/sample_counter
Vout=Vout/sample_counter
Vd=Vd/sample_counter
print("%5i,%10.4f,%10.4f,%10.4f,%10.3f" % (ch,Iout,Vout,Vd,Vin))
#-------Temperature change
#Experimental temperature 1
#R_SCB1=1000*Vd/Iout# RSCB kohm Iout uA
R_SCB1=1000*Vref/Iout# RSCB kohm Iout uA
Temp1 = Temperature(1000*R_SCB1)
Temp1Roots = Temp1.FindRealRoots()
#Experimental temperature 2
#R_SCB2=(10*Vd/(Vout-Vd))# RSCB Kohm
R_SCB2=(10*Vref/(Vout-Vref))# RSCB Kohm
R_SCB2 = 1000*Vout/Iout-10
Temp2 = Temperature(1000*R_SCB2)
Temp2Roots = Temp2.FindRealRoots()
#Theoretical temperature
R_SCBT=(59.948)/(12-Vin)# RSCB kohm
Tempt = Temperature(1000*R_SCBT)
TemptRoots = Tempt.FindRealRoots()
#--------error calculation in the measurements
#Error experimental temperature 1
R_error_e1=ROOT.Math.sqrt(1+(1E6*(Vd/Iout)*(Vd/Iout)))
R_error_e1=(10/Iout)*R_error_e1
T_error_e1= Temperature(R_error_e1+10000)
T_error_e1Roots= T_error_e1.FindRealRoots()
#Error experimental temperature 2
A = 4E-2*((Vd*Vd)/((Vout-Vd)*(Vout-Vd)))
B = 10E-10*((Vout*Vout*10000*10000)/((Vout-Vd)*(Vout-Vd)))
C = 1E-8*((Vd*Vd*10000*10000)/((Vout-Vd)*(Vout-Vd)))
R_error_e2=ROOT.Math.sqrt(A+B+C)
T_error_e2= Temperature(R_error_e2+10000)
T_error_e2Roots= T_error_e2.FindRealRoots()
#--------set data multigraphs
if not ch in data1:
gr=ROOT.TGraph()
gr.SetMarkerStyle(20)
data1[ch]=gr
pass
if not ch in data2:
gr=ROOT.TGraph()
gr.SetMarkerStyle(20)
data2[ch]=gr
pass
if not ch in data3:
gr=ROOT.TGraph()
gr.SetMarkerStyle(20)
data3[ch]=gr
pass
if not ch in data4:
gr=ROOT.TGraph()
gr.SetMarkerStyle(20)
data4[ch]=gr
pass
if not ch in data7:
gr=ROOT.TGraph()
gr.SetMarkerStyle(29)
gr.SetMarkerColor(color_map[MuxMap.index(ch)])
legend.AddEntry(gr,"channel: %5i"%(ch))
data7[ch]=gr
pass
#----------save data from -40 to 25 degrees celsius
if (TemptRoots[0] > -45) and (TemptRoots[0] < 25):
#sample registred
points+=1
#output response
data1[ch].AddPoint(Vout,Temp1Roots[0]) # Temperature (VRB) VS Voltage (SCB)
data2[ch].AddPoint(Vout,Temp2Roots[0]) # Temperature (SCB) VS Voltage (SCB)
#error response
data3[ch].AddPoint(TemptRoots[0],T_error_e1Roots[0]) # Temperature Error (VRB) VS Temperature (THEORETICAL)
data4[ch].AddPoint(TemptRoots[0],T_error_e2Roots[0]) # Temperature Error (SCB) VS Temperature (THEORETICAL)
data7[ch].AddPoint(Temp1Roots[0],Temp1Roots[0]-Temp2Roots[0])##delta Temperature (VRB-SCB) VS Temperature (VRB)
h9.Fill(Temp1Roots[0]-Temp2Roots[0]) #histogram
h10.Fill(Vd)
#storage maximum error on the measure and the maximum found
if channel_error1[ch] < T_error_e1Roots[0]:
channel_error1[ch] = T_error_e1Roots[0]
if channel_error2[ch] < T_error_e1Roots[0]:
channel_error2[ch] = T_error_e2Roots[0]
if maximum_error < (abs(Temp1Roots[0]-Temp2Roots[0])):
maximum_error=abs(Temp1Roots[0]-Temp2Roots[0])
#-------------reset variables for the next cycle
Iout=0
Vout=0
Vd=0
sample_counter=0
flag_data=False
#-------------change in Vin from minmum limit to maximum limit and vice versa
if current_Vin != Vin:
flag_data=True
current_Vin = Vin
pass
#---------------------------------------------------FILL AND FIT DATA--------------------------------------------------
for ch in data1:
mg1.Add(data1[ch])
tf1[ch] = ROOT.TF1('tf1_1_%i'%ch,'pol2')
data1[ch].Fit(tf1[ch],"","",1,2.5)
pass
for ch in data2:
mg2.Add(data2[ch])
tf2[ch] = ROOT.TF1('tf1_2_%i'%ch,'pol2')
data2[ch].Fit(tf2[ch],"","",1,2.5)
pass
for ch in data3:
mg3.Add(data3[ch])
pass
for ch in data4:
mg4.Add(data4[ch])
pass
for ch in data7:
mg5.Add(data7[ch],"PL")
pass
#--------------------------------------------PLOT AND EXPORT DATA------------------------------------------------------------
#------------------DATA 1,2,3,4------------------
#draw the Graph
c1=ROOT.TCanvas("c1","c1",1600,1600)
c1.Divide(2,2)
#first path
c1.cd(1)
mg1.Draw("AP")
mg1.GetXaxis().SetRangeUser(1,2.5)
mg1.GetYaxis().SetRangeUser(-80,60)
#second path
c1.cd(2)
mg2.Draw("AP")
mg2.GetXaxis().SetRangeUser(1,2.5)
mg2.GetYaxis().SetRangeUser(-80,60)
#third path
c1.cd(3)
mg3.Draw("AP")
mg3.GetYaxis().SetRangeUser(0,0.1)
mg3.GetXaxis().SetRangeUser(-80,60)
#fourth path
c1.cd(4)
mg4.Draw("AP")
mg4.GetYaxis().SetRangeUser(0,0.1)
mg4.GetXaxis().SetRangeUser(-80,60)
#print the draw
c1.Update()
c1.SaveAs("C1.pdf")
#--------------------------HISTOGRAM FOR TRANSFER FUNCTION PARAMETERS
for ch in sorted(tf1):
h1.Fill(tf1[ch].GetParameter(0))
h2.Fill(tf1[ch].GetParameter(1))
h3.Fill(tf1[ch].GetParameter(2))
h4.Fill(tf2[ch].GetParameter(0))
h5.Fill(tf2[ch].GetParameter(1))
h6.Fill(tf2[ch].GetParameter(2))
pass
c2=ROOT.TCanvas("c2","c2",1600,1600)
c2.Divide(3,2)
c2.cd(1)
h1.Draw()
c2.cd(2)
h2.Draw()
c2.cd(3)
h3.Draw()
c2.cd(4)
h4.Draw()
c2.cd(5)
h5.Draw()
c2.cd(6)
h6.Draw()
c2.Update()
c2.SaveAs("c2.pdf")
#--------------------------------------DATA 5 AND 6------------------
#--------Canvas file 3
c3=ROOT.TCanvas("c1","c1",1600,800)
c3.Divide(3,1)
#---------error in every channel
c3.cd(1)
mg5.Draw("A pcm plc")#"AP"
legend.Draw()
c3.SetGrid(10,10)
mg5.GetYaxis().SetRangeUser(-0.4,0.4)
mg5.GetXaxis().SetRangeUser(-40,25)
#--------error histogram
c3.cd(2)
h9.Draw()
#---------Vref histogram
c3.cd(3)
h10.Draw()
c3.Update()
c3.SaveAs("c3.pdf")
#-----------------------------------------QUANTILES-------------------------------------------
quantiles=[0.05,0.25,0.50,0.75,0.95]
for rep in range(4):
prob=np.array(quantiles[rep])
y=0.
q=np.array([0.])
y=h9.GetQuantiles(1,q,prob)
quantiles[rep]=q
pass
#-----save_file
wf=ROOT.TFile('out.root','RECREATE')
wf.cd()
mg1.Write()
mg2.Write()
mg3.Write()
mg4.Write()
mg5.Write()
h1.Write()
h2.Write()
h3.Write()
h4.Write()
h5.Write()
h6.Write()
h9.Write()
h10.Write()
wf.Close()
#---------------------------------------Final report-----------------------------------------
#write the output transfer function
fw=open("Transfer_function_report.txt","w")
ss = "ID,DATE,EXAMINATOR,MEAN (ERROR),STDV(ERROR),MAX ERROR,# POINTS, Q(5%), Q(25%), Q(50%), Q(75%), Q(95%), A(MEAN),A(STDV),B(MEAN),B(STDV),C(MEAN),C(STDV),Vref(MEAN),VREF(STDV)\n"
fw.write(ss)
#---data characteristics--------------------------
ss =args.SCB # ID
ss+=","
today = datetime.today()
formatted_date = today.strftime("%Y %m %d")
ss+=formatted_date #DATE
ss+=","
ss+=args.examinator #EXAMINATOR
ss+=","
#----error analysis-------------------------------
ss+="%10.4f,%10.4f,%10.4f,%10i," % (h9.GetMean(1),h9.GetStdDev(1),maximum_error,points)
ss+="%10.4f,%10.4f,%10.4f,%10.4f,%10.4f," % (quantiles[0],quantiles[1],quantiles[2],quantiles[3],quantiles[4])
#----transfer function coefficients---------------
ss+="%10.4f,%10.4f," % (h1.GetMean(1),h1.GetStdDev(1)) # a mean and a stdv
ss+="%10.4f,%10.4f," % (h2.GetMean(1),h2.GetStdDev(1))# b mean and b stdv
ss+="%10.4f,%10.4f," % (h3.GetMean(1),h3.GetStdDev(1))# c mean and c stdv
ss+="%10.4f,%10.20f\n" % (h10.GetMean(1),h10.GetStdDev(1))# Vref mean and Vref stdv
fw.write(ss)
fw.close()
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment