From 43c506ae0cce92ab80add105f8e319d9d7c5e1c3 Mon Sep 17 00:00:00 2001
From: fuchenxi <1256257282@qq.com>
Date: Sun, 22 Sep 2024 17:03:03 +0800
Subject: [PATCH] =?UTF-8?q?=E7=94=B5=E5=AD=90=E5=AD=A6=E5=A4=A7=E9=87=8D?=
 =?UTF-8?q?=E6=9E=84?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 elec/ele_readout.py | 296 +++++++++++++++++---------------------------
 1 file changed, 115 insertions(+), 181 deletions(-)
 mode change 100755 => 100644 elec/ele_readout.py

diff --git a/elec/ele_readout.py b/elec/ele_readout.py
old mode 100755
new mode 100644
index b08e1f3..fdf62f5
--- a/elec/ele_readout.py
+++ b/elec/ele_readout.py
@@ -4,9 +4,9 @@
 '''
 Description: 
     Simulate induced current through BB or CSA amplifier 
-@Date       : 2021/09/02 14:11:57
-@Author     : tanyuhang
-@version    : 1.0
+@Date       : 2024/09/22 15:24:33
+@Author     : tanyuhang, Chenxi Fu
+@version    : 2.0
 '''
 
 import math
@@ -15,10 +15,10 @@ import json
 import ROOT
 
 from current.cal_current import CalCurrent
+from util.math import signal_convolution
 
-# TODO: rewriting in progress
+time_step = 50e-12
 
-# CSA and BB amplifier simulation
 class Amplifier:
     """Get current after amplifier with convolution, for each reading electrode
 
@@ -26,9 +26,11 @@ class Amplifier:
     ---------
     my_current : CalCurrent
         The object of CalCurrent, with induced current and time information
+
     amplifier_name : str
-        The name of amplifier, CSA or BB
-    mintstep : float
+        The name of amplifier
+
+    time_step : float
         The readout time step (bin width)
 
     Attributes
@@ -41,31 +43,28 @@ class Amplifier:
     amplifier_define
         Define parameters and the responce function of amplifier
 
-    sampling_charge
-        Sampling the induced current with readout time step
-
     amplifier_simulation
-        Convolute the induced current with the responce function of amplifier
+        Convolute the induced current with the responce function of amplifier, then scale
 
     Last Modified
     ---------
         2024/09/14
     """
-    def __init__(self, my_current: CalCurrent, amplifier_name: str, mintstep="50e-12"):
+    def __init__(self, my_current: CalCurrent, amplifier_name: str, time_step = time_step):
         self.ele = []
 
         ele_json = "./setting/electronics/" + amplifier_name + ".json"
         with open(ele_json) as f:
-            amplifier_parameters = json.load(f)
+            self.amplifier_parameters = json.load(f)
 
-        self.ele_name = amplifier_parameters['ele_name']
+        self.ele_name = self.amplifier_parameters['ele_name']
         self.read_ele_num = my_current.read_ele_num
-            # each reading electrode has an induced current
-        self.amplifier_define(amplifier_parameters)
-        self.sampling_charge(my_current, mintstep)
-        self.amplifier_simulation()
 
-    def amplifier_define(self, amplifier_parameters: dict):
+        self.amplifier_define()
+        self.fill_amplifier_output(my_current, time_step)
+        self.set_scope_output(my_current)
+
+    def amplifier_define(self):
         """
         Description:
             The parameters of CSA and BB amplifier.
@@ -74,170 +73,105 @@ class Amplifier:
         ---------
             2021/09/09
         """
-        self.CDet = amplifier_parameters['CDet']
+        CDet = self.amplifier_parameters['CDet']
 
-        if amplifier_parameters['ele_name'] == 'CSA':
+        if self.amplifier_parameters['ele_name'] == 'CSA':
             """ CSA parameter initialization"""
-            self.t_rise    = amplifier_parameters['t_rise']
-            self.t_fall    = amplifier_parameters['t_fall']
-            self.trans_imp = amplifier_parameters['trans_imp']
-
-            t_rise = self.t_rise
-            t_fall = self.t_fall
-            self.tau_rise = t_rise/2.2*1e-9
-            self.tau_fall = t_fall/2.2*1e-9
-            if (self.tau_rise == self.tau_fall):
-                self.tau_rise *= 0.9
-            self.sh_max = 0.0  
-
-            self.fill_amplifier_output = self.fill_amplifier_output_CSA
-            self.set_scope_output = self.set_scope_output_CSA
-
-        elif amplifier_parameters['ele_name'] == 'BB':
-            """ BB parameter initialization"""
-            self.BBW       = amplifier_parameters['BBW']
-            self.BBGain    = amplifier_parameters['BBGain']
-            self.BB_imp    = amplifier_parameters['BB_imp']
-            self.OscBW     = amplifier_parameters['OscBW'] 
-
-            tau_C50 = 1.0e-12*50.*self.CDet          #Oscil. RC
-            tau_BW = 0.35/(1.0e9*self.OscBW)/2.2      #Oscil. RC
-            tau_BB_RC = 1.0e-12*self.BB_imp*self.CDet     #BB RC
-            tau_BB_BW = 0.35/(1.0e9*self.BBW)/2.2    #BB Tau
-            self.tau_scope = math.sqrt(pow(tau_C50,2)+pow(tau_BW,2))
-            self.tau_BBA = math.sqrt(pow(tau_BB_RC,2)+pow(tau_BB_BW,2))
-
-            self.fill_amplifier_output = self.fill_amplifier_output_BB
-            self.set_scope_output = self.set_scope_output_BB
-
-    def sampling_charge(self, my_current: CalCurrent, mintstep: float):
-        """ Transform current to charge 
-        with changing bin width to oscilloscope bin width
-        """
-        self.max_num=[]
-        self.itot=[]
-        for i in range(self.read_ele_num):
-            self.max_num.append(my_current.sum_cu[i].GetNbinsX())
-            self.itot.append([0.0]*self.max_num[i])
 
-        self.max_hist_num = my_current.n_bin
-        self.undersampling = int(float(mintstep)/my_current.t_bin)
-        self.time_unit = my_current.t_bin*self.undersampling
-        self.CDet_j = 0     # CSA readout mode
-        
-        self.qtot = [0.0]*self.read_ele_num
-        # self.qtot = [0.0]
-        # get total charge
-        for k in range(self.read_ele_num):
-            i=0
-            for j in range(0,self.max_hist_num,self.undersampling):
-                self.itot[k][i] = my_current.sum_cu[k].GetBinContent(j)
-                self.qtot[k] = self.qtot[k] + self.itot[k][i]*self.time_unit
-                i+=1
-
-    def amplifier_simulation(self):
-        """
-        Description:
-            CSA and BB amplifier Simulation         
-        Parameters:
-        ---------
-        arg1 : int
+            mode = 0
+
+            def pulse_responce_CSA(t):
+                if t < 0: # step function
+                    return 0
+
+                t_rise   = self.amplifier_parameters['t_rise']
+                t_fall   = self.amplifier_parameters['t_fall']
+
+                tau_rise = t_rise/2.2*1e-9
+                tau_fall = t_fall/2.2*1e-9
+                if (tau_rise == tau_fall):
+                    tau_rise *= 0.9
+
+                return tau_fall/(tau_fall+tau_rise) * (math.exp(-t/tau_fall)-math.exp(-t/tau_rise))
+
+            def scale_CSA(output_Q_max, input_Q_tot):
+                """ CSA scale function"""
+                trans_imp = self.amplifier_parameters['trans_imp']
+                Ci = 3.5e-11  #fF
+                Qfrac = 1.0/(1.0+self.CDet*1e-12/Ci)
+
+                if output_Q_max == 0.0:
+                    scale = 0.0
             
-        @Modify:
-        ---------
-            2021/09/09
-        """
-        max_hist_num = int(self.max_hist_num/self.undersampling)
-
-        # Variable Initialization
-        if self.ele_name == 'CSA':
-            IintTime = 2.0*(self.t_rise+self.t_fall)*1e-9/self.time_unit
-            IMaxSh = int(max_hist_num + IintTime)
-            self.shaper_out_Q = [0.0]*IMaxSh
-        elif self.ele_name == 'BB':
-            IintTime = 3.0*self.tau_BBA/self.time_unit
-            IMaxSh = int(max_hist_num + IintTime)
-            self.Iout_BB_RC = [0.0]*IMaxSh
-            self.Iout_C50 = [0.0]*IMaxSh   
-            self.BBGraph = [0.0]*IMaxSh
-
-        self.Vout_scope = [0.0]*IMaxSh
-
-        preamp_Q = [] 
-        for i in range(self.read_ele_num):
-            preamp_Q.append([0.0]*IMaxSh)
-
-        # step for convolution delay
-        step = 1
-
-        for k in range(self.read_ele_num):
-            for i in range(IMaxSh-step):
-                if(i>0 and i <self.max_hist_num-step):
-                    preamp_Q[k][i] = 0.0
-                    for il in range(i,i+step):
-                        preamp_Q[k][i] += self.itot[k][il]*self.time_unit
-                elif (i != 0):
-                    preamp_Q[k][i]=0.0
-
-        for k in range(self.read_ele_num):
-            for i in range(IMaxSh-step):
-                if i >= step:
-                    dif_shaper_Q = preamp_Q[k][i]
+                if mode == 0:
+                    scale = trans_imp * 1e15 * input_Q_tot * Qfrac / output_Q_max     
+                    # scale = trans_imp/(self.CDet*1e-12) #C_D=3.7pF   
+                elif mode == 1:
+                    scale = trans_imp * 1e15 * input_Q_tot / output_Q_max
+
+                return scale
+
+            self.pulse_responce = pulse_responce_CSA
+            self.scale = scale_CSA
+
+        elif self.amplifier_parameters['ele_name'] == 'BB':
+            """ BB parameter initialization"""
+
+            mode = "scope"
+
+            def pulse_responce_BB(t):
+                if t < 0: # step function
+                    return 0
+                
+                BBW       = self.amplifier_parameters['BBW']
+                BB_imp    = self.amplifier_parameters['BB_imp']
+                OscBW     = self.amplifier_parameters['OscBW']   
+                
+                if mode == "scope":
+                    tau_C50 = 1.0e-12 * 50. * CDet          #Oscil. RC
+                    tau_BW = 0.35 / (1.0e9*OscBW) / 2.2      #Oscil. RC
+                    tau_scope = math.sqrt(pow(tau_C50,2)+pow(tau_BW,2))
+
+                    return 1/tau_scope * math.exp(-t/tau_scope)
+
+                elif mode == "RC":
+                    tau_BB_RC = 1.0e-12 * BB_imp * CDet     #BB RC
+                    tau_BB_BW = 0.35 / (1.0e9*BBW) / 2.2    #BB Tau
+                    tau_BBA = math.sqrt(pow(tau_BB_RC,2)+pow(tau_BB_BW,2))
+
+                    return 1/tau_BBA * math.exp(-t/tau_BBA)
+                
                 else:
-                    dif_shaper_Q = 0
-                for j in range(IMaxSh-i):
-                    self.fill_amplifier_output(i, j, dif_shaper_Q)
-                self.set_scope_output(i, k)
-            self.fill_th1f(k, IMaxSh)
-
-    def fill_amplifier_output_CSA(self, i, j, dif_shaper_Q: float):
-        """ Fill CSA out signal, charge"""   
-        self.shaper_out_Q[i+j] += self.tau_fall/(self.tau_fall+self.tau_rise) \
-                                  * dif_shaper_Q*(math.exp(-j*self.time_unit
-                                  / self.tau_fall)-math.exp(
-                                  - j*self.time_unit/self.tau_rise))
-
-    def fill_amplifier_output_BB(self, i, j, dif_shaper_Q: float):
-        """ Fill BB out signal, current"""   
-        self.Iout_C50[i+j] += (dif_shaper_Q)/self.tau_scope \
-                                * math.exp(-j*self.time_unit/self.tau_scope)
-        self.Iout_BB_RC[i+j] += (dif_shaper_Q)/self.tau_BBA \
-                                * math.exp(-j*self.time_unit/self.tau_BBA)
-
-    def set_scope_output_CSA(self, i, k):
-        Ci = 3.5e-11  #fF
-        Qfrac = 1.0/(1.0+self.CDet*1e-12/Ci)
-        Q_max = max(self.shaper_out_Q) if max(self.shaper_out_Q) > 0 else min(self.shaper_out_Q)
-        if Q_max == 0.0:
-            self.Vout_scope[i] = 0.0
-        elif self.CDet_j == 0:
-            self.Vout_scope[i] = self.shaper_out_Q[i]*self.trans_imp\
-                                    * 1e15*self.qtot[k]*Qfrac/Q_max     
-            # self.Vout_scope[i] = self.shaper_out_Q[i]*self.trans_imp/(self.CDet*1e-12) #C_D=3.7pF   
-        elif self.CDet_j == 1:
-            self.Vout_scope[i] = self.shaper_out_Q[i]*self.trans_imp\
-                                    * 1e15*self.qtot[k]/Q_max
-            
-    def set_scope_output_BB(self, i, k):
-        self.BBGraph[i] = 1e3 * self.BBGain * self.Iout_BB_RC[i]
-        R_in = 50 # the input impedance of the amplifier
-        self.Vout_scope[i] = R_in * self.Iout_C50[i]
- 
-    def fill_th1f(self, k, IMaxSh):
-        """ Change amplifier outputs
-            to oscilloscope amplitude [mV]
-            and save in the TH1F
-        """
-        self.ele.append(ROOT.TH1F("electronics %s"%(self.ele_name)+str(k+1), "electronics %s"%(self.ele_name),
-                                IMaxSh, 0, IMaxSh*self.time_unit))
-        # get the max absolute value of the shaper output
-        
-        for i in range(IMaxSh):
-            self.ele[k].SetBinContent(i,self.Vout_scope[i])
-        #Print the max current time of CSA
-        V_max = max(self.Vout_scope) if max(self.Vout_scope) > 0 else min(self.Vout_scope)
-        t_max = self.Vout_scope.index(V_max)
-        print("peak time={:.2e}".format(t_max*self.time_unit))
-
-    def __del__(self):
-        pass
+                    raise NameError(mode,"mode is not defined")
+                
+            def scale_BB(output_Q_max, input_Q_tot):
+                """ BB scale function"""
+
+                if mode == "scope":
+                    R_in = 50
+                    return R_in
+
+                elif mode == "RC":
+                    BBGain = self.amplifier_parameters['BBGain']
+                    return BBGain * 1e3
+                
+            self.pulse_responce = pulse_responce_BB
+            self.scale = scale_BB
+
+    def fill_amplifier_output(self, my_current, time_step):
+        for i in range(self.read_ele_num):
+            sum_cu = my_current.sum_cu[i]
+            n_bin = sum_cu.GetNbinsX()
+            t_bin = sum_cu.GetBinWidth(0)
+            time_duration = n_bin * t_bin
+            self.ele.append(ROOT.TH1F("electronics %s"%(self.ele_name)+str(i+1), "electronics %s"%(self.ele_name),
+                                int(time_duration/time_step), 0, time_duration))
+            self.ele[i].Reset()
+            signal_convolution(sum_cu, self.pulse_responce, self.ele[i])
+    
+    def set_scope_output(self, my_current):
+        for i in range(self.read_ele_num):
+            sum_cu = my_current.sum_cu[i]
+            input_Q_tot = sum_cu.Integral()
+            output_Q_max = self.ele[i].GetMaximum()
+            self.ele[i].Scale(self.scale(output_Q_max, input_Q_tot))
-- 
GitLab