__main__.py

@@ -59,11 +59,11 @@ parser_gen_signal.add_argument('-amp', '--amplifier', type=str, help='amplifier'
 parser_gen_signal.add_argument('-s', '--scan', type=int, help='instance number for scan mode')
 parser_gen_signal.add_argument('--job', type=int, help='flag of run in job')
 
-parser_spaceres = subparsers.add_parser('spaceres', help='space resolution calculation')
-parser_spaceres.add_argument('label', help='LABEL to identify spaceres files')
+parser_telescope = subparsers.add_parser('telescope', help='telescope')
+parser_telescope.add_argument('label', help='LABEL to identify telescope files')
 
-parser_spaceres = subparsers.add_parser('timeres', help='time resolution calculation')
-parser_spaceres.add_argument('det_name', help='name of the detector')
+parser_telescope = subparsers.add_parser('resolution', help='resolution calculation for time, space and energy')
+parser_telescope.add_argument('det_name', help='name of the detector')
 
 parser_tct = subparsers.add_parser('tct', help='TCT simulation')
 parser_tct.add_argument('label', help='LABEL to identify TCT options')
@@ -72,6 +72,7 @@ parser_tct.add_argument('laser', help='name of the laser')
 parser_tct.add_argument('-vol', '--voltage', type=str, help='bias voltage')
 parser_tct.add_argument('-amp', '--amplifier', type=str, help='amplifier')
 parser_tct.add_argument('-s', '--scan', type=int, help='instance number for scan mode')
+parser_tct.add_argument('--job', type=int, help='flag of run in job')
 
 parser_bmos = subparsers.add_parser('bmos', help='Beam Monitor Online System')
 parser_bmos.add_argument('label', help='LABEL to identify BMOS simulations')
@@ -86,14 +87,10 @@ if len(sys.argv) == 1:
 
 kwargs = vars(args)
 
-submodules = ['asic', 'cce', 'cflm', 'current', 'draw', 'elec', 'field', 'fpga', 'gen_signal', 'particle', 'spaceres', 'tct', 'timeres', 'bmos']
-
 submodule = kwargs['subparser_name']
-if submodule not in submodules:
-    raise NameError(submodule)
 
 if kwargs['batch'] != 0:
-    batch_number = kwargs['batch']
+    batch_level = kwargs['batch']
     import re
     from util import batchjob
     destination = submodule
@@ -101,7 +98,8 @@ if kwargs['batch'] != 0:
     command = command.replace('--batch ', '')
     for bs in re.findall('-b* ', command):
         command = command.replace(bs, '')
-    batchjob.main(destination, command, batch_number, args)
+    is_test = vars(args)['test'] 
+    batchjob.main(destination, command, batch_level, is_test)
 else:
     submodule = importlib.import_module(submodule)
     submodule.main(kwargs)

current/cal_current.py

@@ -19,6 +19,7 @@ ROOT.gROOT.SetBatch(True)
 
 from .model import Material
 from particle.carrier_list import CarrierListFromG4P
+from particle.carrier_list import StripCarrierListFromG4P
 from util.math import Vector, signal_convolution
 from util.output import output
 import csv
@@ -293,25 +294,81 @@ class CalCurrent:
                 test_n.Reset()
 
 
-    def save_current(self, my_d, key=None):
-        
-        path = output(__file__, my_d.det_name)
-        if key==None:
-            key = ""
-
-        time = array('d', [999.])
-        current = array('d', [999.])
-        fout = ROOT.TFile(os.path.join(path, "sim-current"+str(key))  + ".root", "RECREATE")
-        t_out = ROOT.TTree("tree", "signal")
-        t_out.Branch("time", time, "time/D")
-        for i in range(self.read_ele_num):
-            t_out.Branch("current"+str(i), current, "current"+str(i)+"/D")
-            for j in range(self.n_bin):
-                current[0]=self.sum_cu[i].GetBinContent(j)
-                time[0]=j*self.t_bin
-                t_out.Fill()
-        t_out.Write()
-        fout.Close()
+    def draw_currents(self, path, tag=""):
+        """
+        @description:
+            Save current in root file
+        @param:
+            None     
+        @Returns:
+            None
+        @Modify:
+            2021/08/31
+        """
+        for read_ele_num in range(self.read_ele_num):
+            c=ROOT.TCanvas("c","canvas1",1600,1300)
+            c.cd()
+            c.Update()
+            c.SetLeftMargin(0.25)
+            # c.SetTopMargin(0.12)
+            c.SetRightMargin(0.15)
+            c.SetBottomMargin(0.17)
+            ROOT.gStyle.SetOptStat(ROOT.kFALSE)
+            ROOT.gStyle.SetOptStat(0)
+
+            #self.sum_cu.GetXaxis().SetTitleOffset(1.2)
+            #self.sum_cu.GetXaxis().SetTitleSize(0.05)
+            #self.sum_cu.GetXaxis().SetLabelSize(0.04)
+            self.sum_cu[read_ele_num].GetXaxis().SetNdivisions(510)
+            #self.sum_cu.GetYaxis().SetTitleOffset(1.1)
+            #self.sum_cu.GetYaxis().SetTitleSize(0.05)
+            #self.sum_cu.GetYaxis().SetLabelSize(0.04)
+            self.sum_cu[read_ele_num].GetYaxis().SetNdivisions(505)
+            #self.sum_cu.GetXaxis().CenterTitle()
+            #self.sum_cu.GetYaxis().CenterTitle() 
+            self.sum_cu[read_ele_num].GetXaxis().SetTitle("Time [s]")
+            self.sum_cu[read_ele_num].GetYaxis().SetTitle("Current [A]")
+            self.sum_cu[read_ele_num].GetXaxis().SetLabelSize(0.08)
+            self.sum_cu[read_ele_num].GetXaxis().SetTitleSize(0.08)
+            self.sum_cu[read_ele_num].GetYaxis().SetLabelSize(0.08)
+            self.sum_cu[read_ele_num].GetYaxis().SetTitleSize(0.08)
+            self.sum_cu[read_ele_num].GetYaxis().SetTitleOffset(1.2)
+            self.sum_cu[read_ele_num].SetTitle("")
+            self.sum_cu[read_ele_num].SetNdivisions(5)
+            self.sum_cu[read_ele_num].Draw("HIST")
+            self.positive_cu[read_ele_num].Draw("SAME HIST")
+            self.negative_cu[read_ele_num].Draw("SAME HIST")
+            self.gain_positive_cu[read_ele_num].Draw("SAME HIST")
+            self.gain_negative_cu[read_ele_num].Draw("SAME HIST")
+            self.sum_cu[read_ele_num].Draw("SAME HIST")
+
+            self.positive_cu[read_ele_num].SetLineColor(877)#kViolet-3
+            self.negative_cu[read_ele_num].SetLineColor(600)#kBlue
+            self.gain_positive_cu[read_ele_num].SetLineColor(617)#kMagneta+1
+            self.gain_negative_cu[read_ele_num].SetLineColor(867)#kAzure+7
+            self.sum_cu[read_ele_num].SetLineColor(418)#kGreen+2
+
+            self.positive_cu[read_ele_num].SetLineWidth(2)
+            self.negative_cu[read_ele_num].SetLineWidth(2)
+            self.gain_positive_cu[read_ele_num].SetLineWidth(2)
+            self.gain_negative_cu[read_ele_num].SetLineWidth(2)
+            self.sum_cu[read_ele_num].SetLineWidth(2)
+            c.Update()
+
+            legend = ROOT.TLegend(0.5, 0.2, 0.8, 0.5)
+            legend.AddEntry(self.negative_cu[read_ele_num], "electron", "l")
+            legend.AddEntry(self.positive_cu[read_ele_num], "hole", "l")
+            legend.AddEntry(self.gain_negative_cu[read_ele_num], "gain electron", "l")
+            legend.AddEntry(self.gain_positive_cu[read_ele_num], "gain hole", "l")
+            legend.AddEntry(self.sum_cu[read_ele_num], "e+h", "l")
+            legend.SetBorderSize(0)
+            #legend.SetTextFont(43)
+            legend.SetTextSize(0.08)
+            legend.Draw("same")
+            c.Update()
+            c.SaveAs(path+'/'+tag+"No_"+str(read_ele_num+1)+"electrode"+"_basic_infor.pdf")
+            c.SaveAs(path+'/'+tag+"No_"+str(read_ele_num+1)+"electrode"+"_basic_infor.root")
+            del c
 
     
 class CalCurrentGain(CalCurrent):
@@ -494,64 +551,3 @@ class CalCurrentLaser(CalCurrent):
             self.gain_positive_cu[i] = convolved_gain_positive_cu
             self.gain_negative_cu[i] = convolved_gain_negative_cu
             self.sum_cu[i] = convolved_sum_cu
-
-class StripCarrierListFromG4P:
-    def __init__(self, material, my_g4p, batch):
-        if (material == "SiC"):
-            self.energy_loss = 8.4 #ev
-        elif (material == "Si"):
-            self.energy_loss = 3.6 #ev
-
-        if batch == 0:
-            h1 = ROOT.TH1F("Edep_device", "Energy deposition in Detector", 100, 0, max(my_g4p.edep_devices)*1.1)
-            for i in range (len(my_g4p.edep_devices)):
-                h1.Fill(my_g4p.edep_devices[i])
-            max_event_bin=h1.GetMaximumBin()
-            bin_wide=max(my_g4p.edep_devices)*1.1/100
-            for j in range (len(my_g4p.edep_devices)):
-                #compare to experimental data
-                if (my_g4p.edep_devices[j]<0.084 and my_g4p.edep_devices[j]>0.083):
-                    try_p=1
-                    for single_step in my_g4p.p_steps_current[j]:
-                        if abs(single_step[0]-my_g4p.p_steps_current[j][0][0])>5:
-                            try_p=0
-                    if try_p==1:
-                        self.batch_def(my_g4p,j)
-                        batch = 1
-                        break
-
-            if batch == 0:
-                print("the sensor didn't have particles hitted")
-                raise ValueError
-        else:
-            self.batch_def(my_g4p,batch)
-
-    def batch_def(self,my_g4p,j):
-        self.beam_number = j
-        self.track_position = [[single_step[0],single_step[1],single_step[2],1e-9] for single_step in my_g4p.p_steps_current[j]]
-        self.tracks_step = my_g4p.energy_steps[j]
-        self.tracks_t_energy_deposition = my_g4p.edep_devices[j] #为什么不使用？
-        self.ionized_pairs = [step*1e6/self.energy_loss for step in self.tracks_step]
-
-# TODO: change this to a method of CalCurrent
-def save_current(my_d,my_l,my_current,my_f,key):
-    if "planar3D" in my_d.det_model or "planarRing" in my_d.det_model:
-        path = os.path.join('output', 'pintct', my_d.det_name, )
-    elif "lgad3D" in my_d.det_model:
-        path = os.path.join('output', 'lgadtct', my_d.det_name, )
-    create_path(path) 
-    L = eval("my_l.{}".format(key))
-    #L is defined by different keys
-    time = array('d', [999.])
-    current = array('d', [999.])
-    fout = ROOT.TFile(os.path.join(path, "sim-TCT-current") + str(L) + ".root", "RECREATE")
-    t_out = ROOT.TTree("tree", "signal")
-    t_out.Branch("time", time, "time/D")
-    for i in range(my_f.read_ele_num):
-        t_out.Branch("current"+str(i), current, "current"+str(i)+"/D")
-        for j in range(my_current.n_bin):
-            current[0]=my_current.sum_cu[i].GetBinContent(j)
-            time[0]=j*my_current.t_bin
-            t_out.Fill()
-        t_out.Write()
-        fout.Close()
\ No newline at end of file

current/cal_current_diffuse.py

@@ -17,6 +17,7 @@ import ROOT
 ROOT.gROOT.SetBatch(True)
 
 from .model import Material
+from particle.carrier_list import PixelCarrierListFromG4P
 from util.math import Vector
 
 t_bin = 50e-12
@@ -229,63 +230,3 @@ class CalCurrentPixel:
                     Hit["charge"].append(charge)       
         return Hit["index"],Hit["charge"]
 
-
-class PixelCarrierListFromG4P:
-    def __init__(self, my_d,my_g4p):
-        """
-        Description:
-            Events position and energy depositon
-        Parameters:
-            material : string
-                deciding the energy loss of MIP
-            my_g4p : Particles
-            batch : int
-                batch = 0: Single event, select particle with long enough track
-                batch != 0: Multi event, assign particle with batch number
-        Modify:
-            2022/10/25
-        """
-        batch = len(my_g4p.localposition)
-        layer = len(my_d.lt_z)
-        material = my_d.material
-        self.pixelsize_x = my_d.p_x
-        self.pixelsize_y = my_d.p_y
-        self.pixelsize_z = my_d.p_z
-        
-        if (material == "SiC"):
-            self.energy_loss = 8.4 #ev
-        elif (material == "Si"):
-            self.energy_loss = 3.6 #ev
-        
-        self.track_position, self.ionized_pairs= [],[]
-        self.layer= layer
-        for j in range(batch):
-            self.single_event(my_g4p,j)
-
-    def single_event(self,my_g4p,j):
-        s_track_position,s_energy= [],[]
-        for i in range(self.layer):
-            position = []
-            energy = []
-            name = "Layer_"+str(i)
-            #print(name)
-            for k in range(len(my_g4p.devicenames[j])):
-                px,py,pz = self.split_name(my_g4p.devicenames[j][k])
-                if name in my_g4p.devicenames[j][k]:
-                    tp = [0 for i in range(3)]
-                    tp[0] = my_g4p.localposition[j][k][0]+(px-0.5)*self.pixelsize_x
-                    tp[1] = my_g4p.localposition[j][k][1]+(py-0.5)*self.pixelsize_y
-                    tp[2] = my_g4p.localposition[j][k][2]+self.pixelsize_z/2
-                    position.append(tp) 
-                    energy.append(my_g4p.energy_steps[j][k])
-            s_track_position.append(position)
-            pairs = [step*1e6/self.energy_loss for step in energy]
-            s_energy.append(pairs)
-            del position,energy
-        self.track_position.append(s_track_position)
-        self.ionized_pairs.append(s_energy)
-        
-    def split_name(self,volume_name):
-        parts = volume_name.split('_')
-        return int(parts[1]),int(parts[2]),int(parts[4])
-

elec/__init__.py

@@ -4,37 +4,14 @@ from util.output import output
 
 def main(kwargs):
     label = kwargs['label'] # Operation label or detector name
-    name = kwargs['name']
-    tct = kwargs['tct']
+    name = kwargs['name'] # readout electronics name
     os.makedirs('output/elec/{}'.format(name), exist_ok=True)
 
     if label == 'trans':
         subprocess.run(['ngspice -b setting/electronics/{}.cir'.format(name)], shell=True)
-    elif label == 'get_fig':
-        from . import ngspice_get_fig
-        file_path = output(__file__)
-        ngspice_get_fig.main(label, file_path)
     elif label == 'readout':
         from . import readout
-        readout.main(label)
+        readout.main(name)
     else:
-        if tct == None:
-            from . import ngspice_set_input
-            from . import ngspice_set_tmp_cir
-            input_p = ngspice_set_input.set_input(label)
-            input_c=','.join(input_p)
-            ngspice_set_tmp_cir.ngspice_set_tmp_cir(input_c, label, name)
-            subprocess.run(['ngspice -b output/elec/{}/{}_tmp.cir'.format(label, name)], shell=True)
-            file_path = output(__file__, label)
-            from . import ngspice_get_fig
-            ngspice_get_fig.main(name, file_path)
-        if tct != None:
-            from . import ngspice_set_input
-            from . import ngspice_set_tmp_cir
-            input_p = ngspice_set_input.set_input(label, tct)
-            input_c=','.join(input_p)
-            ngspice_set_tmp_cir.ngspice_set_tmp_cir(input_c, label, name, tct)
-            subprocess.run(['ngspice -b output/elec/{}/{}{}_tmp.cir'.format(label, name, tct)], shell=True)
-            file_path = output(__file__, label)
-            from . import ngspice_get_fig
-            ngspice_get_fig.main(name, file_path, tct)
+        raise NameError
+

elec/ngspice_get_fig.py

-#!/usr/bin/env python3
-
-import sys
-import os
-
-import numpy
-import ROOT
-ROOT.gROOT.SetBatch(True)
-
-from util.output import output
-
-def read_file(file_path,file_name):
-    with open(file_path + '/' + file_name) as f:
-        lines = f.readlines()
-        time,volt = [],[]
-
-        for line in lines:
-            time.append(float(line.split()[0])*1e9)
-            volt.append(float(line.split()[1])*1e3)
-
-    time = numpy.array(time,dtype='float64')
-    volt = numpy.array(volt,dtype='float64')
-
-    return time,volt
-
-def main(elec_name, file_path, key=None):
-    if key is None:
-        key = ''
-    fig_name = os.path.join(file_path, elec_name+key+'.pdf')
-    time,volt = [],[]
-
-    time,volt = read_file(file_path, elec_name+key+'.raw')
-    length = len(time)
-    t_min, t_max = time[0], time[-1]
-
-    ROOT.gROOT.SetBatch()    
-    c = ROOT.TCanvas('c','c',700,600)
-    c.SetMargin(0.2,0.1,0.2,0.1)
-    f1 = ROOT.TGraph(length,time,volt)
-    f1.SetTitle(' ')
-
-    f1.SetLineColor(2)
-    f1.SetLineWidth(2)
-
-    f1.GetXaxis().SetTitle('Time [ns]')
-    f1.GetXaxis().SetLimits(t_min, t_max)
-    f1.GetXaxis().CenterTitle()
-    f1.GetXaxis().SetTitleSize(0.08)
-    f1.GetXaxis().SetLabelSize(0.08)
-    f1.GetXaxis().SetNdivisions(5)
-    f1.GetXaxis().SetTitleOffset(1)
-
-    f1.GetYaxis().SetTitle('Voltage [mV]')
-    f1.GetYaxis().CenterTitle()
-    f1.GetYaxis().SetTitleSize(0.08)
-    f1.GetYaxis().SetLabelSize(0.08)
-    f1.GetYaxis().SetNdivisions(5)
-    f1.GetYaxis().SetTitleOffset(1)
-
-    c.cd()
-    f1.Draw('AL')
-    c.SaveAs(fig_name)
-
-if __name__ == '__main__':
-    main()
-    

elec/ngspice_set_input.py

-import os
-from array import array
-
-import ROOT
-ROOT.gROOT.SetBatch(True)
-
-from util.output import output
-
-def set_input(det_name, key=None):
-    current=[]
-    time=[]
-    if key == None:
-        key = ""
-    path = "output/current/{}".format(det_name)
-    myFile = ROOT.TFile(os.path.join(path, "sim-current"+str(key))+".root")
-
-    myt = myFile.tree
-    for entry in myt:
-       current.append(entry.current0) # current[i], i for electrode number
-       time.append(entry.time)
-    input_c=[]
-    if abs(min(current))>max(current): #set input signal
-        c_max=min(current)
-        for i in range(0, len(current)):
-            if current[i] < c_max * 0.01:
-                input_c.append(str(0))
-                input_c.append(str(0))
-                input_c.append(str(time[i]))
-                input_c.append(str(0))
-                break
-            else:
-                current[i]=0
-        for j in range(i, len(current)):
-            input_c.append(str(time[j]))
-            input_c.append(str(current[j]))
-            if current[j] > c_max * 0.01:
-                break
-        input_c.append(str(time[j]))
-        input_c.append(str(0))
-        input_c.append(str(time[len(time)-1]))
-        input_c.append(str(0))
-        for k in range(j, len(current)):
-            current[i]=0
-    else:
-        c_max=max(current)
-        for i in range(0, len(current)):
-            current[i]=0
-            if current[i] > c_max * 0.01:
-                input_c.append(str(0))
-                input_c.append(str(0))
-                input_c.append(str(time[i]))
-                input_c.append(str(0))
-                break
-        for j in range(i, len(current)):
-            input_c.append(str(time[j]))
-            input_c.append(str(current[j]))
-            if current[j] < c_max * 0.01:
-                break
-        input_c.append(str(time[j]))
-        input_c.append(str(0))
-        input_c.append(str(time[len(time)-1]))
-        input_c.append(str(0))
-        for k in range(j, len(current)):
-            current[i]=0
-    in_put=array("d",[0.])
-    t=array("d",[0.])
-    out_path = output(__file__, det_name)
-    fout = ROOT.TFile(os.path.join(out_path, "input"+str(key))+".root", "RECREATE")
-    t_out = ROOT.TTree("tree", "signal")
-    t_out.Branch("time", t, "time/D")
-    t_out.Branch("current", in_put, "current/D")
-    n_bin = myt.GetEntries()
-    for m in range(n_bin):
-        in_put[0]=current[m]
-        t[0]=time[m]
-        t_out.Fill()
-    t_out.Write()
-    fout.Close()
-    return input_c

elec/ngspice_set_tmp_cir.py

-import re
-def ngspice_set_tmp_cir(input_c, det_name, ele_name, label=None):
-    if label is None:
-        label = ''
-    with open('./setting/electronics/{}.cir'.format(ele_name), 'r') as f:
-        lines = f.readlines()
-        for i in range(len(lines)):
-            if lines[i].startswith('I1'):
-                # replace pulse by PWL
-                lines[i] = re.sub(r"pulse" + r".*", 'PWL('+str(input_c)+') \n', lines[i], flags=re.IGNORECASE)
-            if lines[i].startswith('wrdata'):
-                # replace output file name & path
-                lines[i] = re.sub(r".*" + r".raw", "wrdata output/elec/{}/{}{}.raw".format(det_name, ele_name, label), lines[i])
-        f.close()
-    with open('output/elec/{}/{}{}_tmp.cir'.format(det_name, ele_name, label), 'w+') as f:
-        f.writelines(lines)
-        f.close()
\ No newline at end of file

field/devsim_draw.py

@@ -287,3 +287,4 @@ def draw2D(x,y,value,title,v,path):
     graph_1d.GetYaxis().SetTitle("Potential")
     canvas1.SaveAs(os.path.join(path, title+"{}_1d.png".format(v)))
     canvas1.SaveAs(os.path.join(path, title+"{}_1d.root".format(v)))
+

gen_signal/__init__.py

@@ -5,12 +5,18 @@ def main(kwargs):
 
     if label == 'signal':
         if scan_number != None:
-            if job_number != None:
-                from . import gen_signal_scan
-                gen_signal_scan.job_main(kwargs)
-            else:
-                from . import gen_signal_scan
-                gen_signal_scan.main(kwargs)
+            from util import batchjob
+            scan_number = kwargs['scan']
+            det_name = kwargs['det_name']
+            for i in range(scan_number):
+                args = ['gen_signal', '--job', str(i), det_name]
+                command = ' '.join(args) 
+                print(command)
+                destination = 'gen_signal'
+                batchjob.main(destination, command, 1, is_test=False)
+        elif job_number != None:
+            from . import gen_signal_scan
+            gen_signal_scan.main(kwargs)
         else:
             from . import gen_signal_main
             gen_signal_main.main(kwargs)
\ No newline at end of file

gen_signal/beam_monitor.py

@@ -16,7 +16,7 @@ def get_beam_number(my_g4p,ele_current):
     create_path(path) 
     number = array('d',[999.])
     hittotal = array('d',[999.])
-    number[0] = int(-ele_current.max_BB_height/18.8)
+    number[0] = int(-ele_current.max_Broad_Band_height/18.8)
     hittotal[0]=my_g4p.hittotal
     fout = ROOT.TFile(path + "beam_monitor.root", "RECREATE")
     t_out = ROOT.TTree("tree", "beam_number")

gen_signal/build_device.py

@@ -31,7 +31,7 @@ class Detector:
         self.l_y = self.device_dict['ly']  
         self.l_z = self.device_dict['lz'] 
         
-        self.voltage = self.device_dict['bias']['voltage'] 
+        self.voltage = float(self.device_dict['bias']['voltage'])
         self.temperature = self.device_dict['temperature']
         self.material = self.device_dict['material']
 

gen_signal/draw_save.py

@@ -209,117 +209,6 @@ def draw_drift_path(my_d,my_g4p,my_f,my_current,path):
     c1.SaveAs(path+'/'+my_d.det_model+"_drift_path.root")
     del c1
 
-def draw_current(my_d, my_current, ele_current, read_ele_num, model, path, tag=""):
-    """
-    @description:
-        Save current in root file
-    @param:
-        None     
-    @Returns:
-        None
-    @Modify:
-        2021/08/31
-    """
-    c=ROOT.TCanvas("c","canvas1",1600,1300)
-    c.cd()
-    c.Update()
-    c.SetLeftMargin(0.25)
-    # c.SetTopMargin(0.12)
-    c.SetRightMargin(0.15)
-    c.SetBottomMargin(0.17)
-    ROOT.gStyle.SetOptStat(ROOT.kFALSE)
-    ROOT.gStyle.SetOptStat(0)
-
-    #my_current.sum_cu.GetXaxis().SetTitleOffset(1.2)
-    #my_current.sum_cu.GetXaxis().SetTitleSize(0.05)
-    #my_current.sum_cu.GetXaxis().SetLabelSize(0.04)
-    my_current.sum_cu[read_ele_num].GetXaxis().SetNdivisions(510)
-    #my_current.sum_cu.GetYaxis().SetTitleOffset(1.1)
-    #my_current.sum_cu.GetYaxis().SetTitleSize(0.05)
-    #my_current.sum_cu.GetYaxis().SetLabelSize(0.04)
-    my_current.sum_cu[read_ele_num].GetYaxis().SetNdivisions(505)
-    #my_current.sum_cu.GetXaxis().CenterTitle()
-    #my_current.sum_cu.GetYaxis().CenterTitle() 
-    my_current.sum_cu[read_ele_num].GetXaxis().SetTitle("Time [s]")
-    my_current.sum_cu[read_ele_num].GetYaxis().SetTitle("Current [A]")
-    my_current.sum_cu[read_ele_num].GetXaxis().SetLabelSize(0.08)
-    my_current.sum_cu[read_ele_num].GetXaxis().SetTitleSize(0.08)
-    my_current.sum_cu[read_ele_num].GetYaxis().SetLabelSize(0.08)
-    my_current.sum_cu[read_ele_num].GetYaxis().SetTitleSize(0.08)
-    my_current.sum_cu[read_ele_num].GetYaxis().SetTitleOffset(1.2)
-    my_current.sum_cu[read_ele_num].SetTitle("")
-    my_current.sum_cu[read_ele_num].SetNdivisions(5)
-    my_current.sum_cu[read_ele_num].Draw("HIST")
-    my_current.positive_cu[read_ele_num].Draw("SAME HIST")
-    my_current.negative_cu[read_ele_num].Draw("SAME HIST")
-    my_current.gain_positive_cu[read_ele_num].Draw("SAME HIST")
-    my_current.gain_negative_cu[read_ele_num].Draw("SAME HIST")
-    my_current.sum_cu[read_ele_num].Draw("SAME HIST")
-
-    my_current.positive_cu[read_ele_num].SetLineColor(877)#kViolet-3
-    my_current.negative_cu[read_ele_num].SetLineColor(600)#kBlue
-    my_current.gain_positive_cu[read_ele_num].SetLineColor(617)#kMagneta+1
-    my_current.gain_negative_cu[read_ele_num].SetLineColor(867)#kAzure+7
-    my_current.sum_cu[read_ele_num].SetLineColor(418)#kGreen+2
-
-    my_current.positive_cu[read_ele_num].SetLineWidth(2)
-    my_current.negative_cu[read_ele_num].SetLineWidth(2)
-    my_current.gain_positive_cu[read_ele_num].SetLineWidth(2)
-    my_current.gain_negative_cu[read_ele_num].SetLineWidth(2)
-    my_current.sum_cu[read_ele_num].SetLineWidth(2)
-    c.Update()
-
-    if abs(ele_current[read_ele_num].GetMinimum()) > abs(ele_current[read_ele_num].GetMaximum()):
-        rightmax = 1.1*ele_current[read_ele_num].GetMinimum()
-    else:
-        rightmax = 1.1*ele_current[read_ele_num].GetMaximum()
-
-    if rightmax == 0.0:
-        n_scale = 0.0
-    elif abs(ele_current[read_ele_num].GetMinimum()) > abs(ele_current[read_ele_num].GetMaximum()):
-        n_scale = ROOT.gPad.GetUymin() / rightmax
-    else:
-        n_scale = ROOT.gPad.GetUymax() / rightmax
-    """ 
-    ele_current[read_ele_num].Scale(n_scale)
-    ele_current[read_ele_num].Draw("SAME HIST")
-    ele_current[read_ele_num].SetLineWidth(2)   
-    ele_current[read_ele_num].SetLineColor(8)
-    ele_current[read_ele_num].SetLineColor(2)
-    c.Update()
-
-    axis = ROOT.TGaxis(ROOT.gPad.GetUxmax(), ROOT.gPad.GetUymin(), 
-                       ROOT.gPad.GetUxmax(), ROOT.gPad.GetUymax(), 
-                       min(0,rightmax), max(0,rightmax), 510, "+L")
-    axis.SetLineColor(2)
-    axis.SetTextColor(2)
-    axis.SetTextSize(0.02)
-    axis.SetTextFont(40)
-    axis.SetLabelColor(2)
-    axis.SetLabelSize(0.035)
-    axis.SetLabelFont(42)
-    axis.SetTitle("Ampl [V]")
-    axis.SetTitleFont(40)
-    axis.SetTitleOffset(1.2)
-    #axis.CenterTitle()
-    axis.Draw("SAME HIST") """
-
-    legend = ROOT.TLegend(0.5, 0.2, 0.8, 0.5)
-    legend.AddEntry(my_current.negative_cu[read_ele_num], "electron", "l")
-    legend.AddEntry(my_current.positive_cu[read_ele_num], "hole", "l")
-    #legend.AddEntry(my_current.gain_negative_cu[read_ele_num], "gain electron", "l")
-    #legend.AddEntry(my_current.gain_positive_cu[read_ele_num], "gain hole", "l")
-    legend.AddEntry(my_current.sum_cu[read_ele_num], "e+h", "l")
-    #legend.AddEntry(ele_current, "electronics", "l")
-    legend.SetBorderSize(0)
-    #legend.SetTextFont(43)
-    legend.SetTextSize(0.08)
-    legend.Draw("same")
-    c.Update()
-    c.SaveAs(path+'/'+model+my_d.det_model+tag+"No_"+str(read_ele_num+1)+"electrode"+"_basic_infor.pdf")
-    c.SaveAs(path+'/'+model+my_d.det_model+tag+"No_"+str(read_ele_num+1)+"electrode"+"_basic_infor.root")
-    del c
-
 def cce(my_current, path):
     charge=array('d')
     x=array('d')
@@ -341,23 +230,3 @@ def cce(my_current, path):
     c1.SaveAs(path+"/cce.pdf")
     c1.SaveAs(path+"/cce.root")
 
-def save_signal_time_resolution(my_d,batch_number,sum_cu,ele_current,my_g4p,start_n):
-    """ Save induced current after amplification"""
-
-    output_path = output(__file__, my_d.det_name, 'batch')
-    for k in range(ele_current.read_ele_num):
-        charge = sum_cu[k].Integral()
-        charge_str = "_charge=%.2f_"%(charge*1e15)  #fc
-        e_dep = "dep=%.5f_"%(my_g4p.edep_devices[batch_number-start_n]) #mv
-        output_file = output_path + "/t_" +str(batch_number)+charge_str+e_dep+"events.csv"
-
-        f1 = open(output_file,"w")
-        f1.write("time[ns], Amplitude [mV] \n")
-        for i in range(ele_current.amplified_current[k].GetNbinsX()):
-            f1.write("%s,%s \n"%(i*ele_current.amplified_current[k].GetBinWidth(0),
-                                    ele_current.amplified_current[k][i]))
-        f1.close()
-
-        print("output_file:%s"%output_file)
-
-    del ele_current

gen_signal/gen_signal_main.py

@@ -23,7 +23,7 @@ from particle import g4_time_resolution as g4t
 from field import devsim_field as devfield
 from current import cal_current as ccrt
 from elec import readout as rdo
-from .draw_save import energy_deposition, draw_drift_path, draw_current, cce
+from .draw_save import energy_deposition, draw_drift_path, cce
 from util.output import output
 
 
@@ -69,20 +69,17 @@ def main(kwargs):
     g4_seed = random.randint(0,1e7)
     my_g4p = g4t.Particles(my_d, absorber, g4_seed)
     my_current = ccrt.CalCurrentG4P(my_d, my_f, my_g4p, 0)
+    ele_current = rdo.Amplifier(my_current.sum_cu, amplifier)
 
     now = time.strftime("%Y_%m%d_%H%M%S")
     path = output(__file__, my_d.det_name, now)
-
     #energy_deposition(my_g4p)   # Draw Geant4 depostion distribution
     draw_drift_path(my_d,my_g4p,my_f,my_current,path)
+    my_current.draw_currents(path) # Draw current
+    ele_current.draw_waveform(my_current.sum_cu, path) # Draw waveform
 
-    my_current.save_current(my_d)
-    if 'ngspice' not in amplifier:
-        ele_current = rdo.Amplifier(my_current.sum_cu, amplifier)
-        for i in range(my_current.read_ele_num):
-            draw_current(my_d, my_current, ele_current.amplified_current, i, ele_current.amplified_current_name, path) # Draw current
-        if 'strip' in my_d.det_model:
-            cce(my_current, path)
+    if 'strip' in my_d.det_model:
+        cce(my_current, path)
     
     del my_f
     end = time.time()

gen_signal/gen_signal_main_SiCstrip.py

+#!/usr/bin/env python3
+# -*- encoding: utf-8 -*-
+'''
+@Description: The main program of Raser induced current simulation,(quickly checkout under different Votage)      
+@Date       : 2025/02/11 
+@Author     : zhulin
+@version    : 2.0
+'''
+import sys
+import os
+import array
+import time
+import subprocess
+import json
+import random
+
+import ROOT
+ROOT.gROOT.SetBatch(True)
+import geant4_pybind as g4b
+
+from . import build_device as bdv
+from particle import g4_time_resolution as g4t
+from field import devsim_field as devfield
+from current import cal_current as ccrt
+from elec import readout as rdo
+from .draw_save import energy_deposition, draw_drift_path, draw_current, cce
+from util.output import output
+
+
+def main(kwargs):
+    """
+    Description:
+        The main program of Raser induced current simulation      
+    Parameters:
+    ---------
+    dset : class
+        Parameters of simulation
+    Function or class:
+        Detector -- Define the basic parameters and mesh structure of the detector
+        DevsimField -- Get the electric field and weighting potential 
+        Particles -- Electron and hole paris distibution
+        CalCurrent -- Drift of e-h pais and induced current
+        Amplifier -- Readout electronics simulation  
+    Modify:
+    ---------
+        2021/09/02
+    """
+    start = time.time()
+
+    det_name = kwargs['det_name']
+    my_d = bdv.Detector(det_name)
+    if kwargs['voltage'] != None:
+        if kwargs['absorber'] != None:
+            absorber = kwargs['absorber']
+        else:
+            absorber = my_d.absorber
+        if kwargs['amplifier'] != None:
+            amplifier = kwargs['amplifier']
+        else:
+            amplifier = my_d.amplifier
+    
+        g4_seed = random.randint(0,1e7)
+        my_g4p = g4t.Particles(my_d, absorber, g4_seed)
+
+        voltage_max = int(kwargs['voltage'])
+        for i in range(1,abs(voltage_max)+1):
+            if voltage_max<0:
+                voltage = -1*int(i)
+            else:
+                voltage = int(i)
+            
+            # my_f = devfield.DevsimField(my_d.device, my_d.dimension, voltage, my_d.read_out_contact, my_d.irradiation_flux)
+            if "strip" in det_name:
+                my_f = devfield.DevsimField(my_d.device, my_d.dimension, voltage, my_d.read_ele_num, my_d.l_z)
+            else: 
+                my_f = devfield.DevsimField(my_d.device, my_d.dimension, voltage, 1, my_d.l_z)
+
+            
+            # my_current = ccrt.CalCurrentG4P(my_d, my_f, my_g4p, 0)
+            if "strip" in det_name:
+                my_current = ccrt.CalCurrentStrip(my_d, my_f, my_g4p, 0)
+            else: 
+                my_current = ccrt.CalCurrentG4P(my_d, my_f, my_g4p, 0)
+
+            now = time.strftime("%Y_%m%d_%H%M%S")
+            path = output(__file__, my_d.det_name, now)
+
+            #energy_deposition(my_g4p)   # Draw Geant4 depostion distribution
+            draw_drift_path(my_d,my_g4p,my_f,my_current,path)
+
+            my_current.save_current(my_d)
+            if 'ngspice' not in amplifier:
+                ele_current = rdo.Amplifier(my_current.sum_cu, amplifier)
+                for i in range(my_current.read_ele_num):
+                    draw_current(my_d, my_current, ele_current.amplified_current, i, ele_current.amplified_current_name, path) # Draw current
+                if 'strip' in my_d.det_model:
+                    cce(my_current, path)
+            
+            del my_f
+            end = time.time()
+            print("total_time:%s"%(end-start))
+    
+    else:
+        if kwargs['absorber'] != None:
+            absorber = kwargs['absorber']
+        else:
+            absorber = my_d.absorber
+        if kwargs['amplifier'] != None:
+            amplifier = kwargs['amplifier']
+        else:
+            amplifier = my_d.amplifier
+    
+        g4_seed = random.randint(0,1e7)
+        my_g4p = g4t.Particles(my_d, absorber, g4_seed)
+
+        voltage_max = int(my_d.voltage)
+        for i in range(500,abs(voltage_max)+1,10):
+            if voltage_max<0:
+                # voltage = -1*int(i)
+                voltage = -1*float(i)
+            else:
+                # voltage = int(i)
+                voltage = float(i)
+            
+            my_f = devfield.DevsimField(my_d.device, my_d.dimension, voltage, my_d.read_out_contact, my_d.irradiation_flux)
+            
+            my_current = ccrt.CalCurrentG4P(my_d, my_f, my_g4p, 0)
+            # if "strip" in det_name:
+            #     my_current = ccrt.CalCurrentStrip(my_d, my_f, my_g4p, 0)
+            # else: 
+            #     my_current = ccrt.CalCurrentG4P(my_d, my_f, my_g4p, 0)
+
+            now = time.strftime("%Y_%m%d_%H%M%S")
+            path = output(__file__, my_d.det_name, now)
+
+            #energy_deposition(my_g4p)   # Draw Geant4 depostion distribution
+            draw_drift_path(my_d,my_g4p,my_f,my_current,path)
+
+            my_current.save_current(my_d)
+            if 'ngspice' not in amplifier:
+                ele_current = rdo.Amplifier(my_current.sum_cu, amplifier)
+                for i in range(my_current.read_ele_num):
+                    draw_current(my_d, my_current, ele_current.amplified_current, i, ele_current.amplified_current_name, path) # Draw current
+                if 'strip' in my_d.det_model:
+                    cce(my_current, path)
+            
+            del my_f
+            end = time.time()
+            print("total_time:%s"%(end-start))
+
+
+if __name__ == '__main__':
+    args = sys.argv[1:]
+    kwargs = {}
+    for arg in args:
+        key, value = arg.split('=')
+        kwargs[key] = value
+    main(kwargs)
+    
\ No newline at end of file

gen_signal/gen_signal_scan.py

@@ -22,11 +22,9 @@ from particle import g4_time_resolution as g4t
 from field import devsim_field as devfield
 from current import cal_current as ccrt
 from elec import readout as rdo
-from . import draw_save
 from util.output import output
 
 
-
 def batch_loop(my_d, my_f, my_g4p, amplifier, g4_seed, total_events, instance_number):
     """
     Description:
@@ -56,12 +54,17 @@ def batch_loop(my_d, my_f, my_g4p, amplifier, g4_seed, total_events, instance_nu
             effective_number += 1
             my_current = ccrt.CalCurrentG4P(my_d, my_f, my_g4p, event-start_n)
             ele_current = rdo.Amplifier(my_current.sum_cu, amplifier)
-            draw_save.save_signal_time_resolution(my_d,event,my_current.sum_cu,ele_current,my_g4p,start_n)
+
+            e_dep = "%.5f"%(my_g4p.edep_devices[event-start_n]) #mv
+            tag = "event" + str(event) + "_" + "Edep" + str(e_dep)
+            # need to add ngspice branch 
+            ele_current.save_signal_TTree(output(__file__, my_d.det_name, 'batch'), tag)
             del ele_current
+
     detection_efficiency =  effective_number/(end_n-start_n) 
     print("detection_efficiency=%s"%detection_efficiency)
 
-def job_main(kwargs):
+def main(kwargs):
     det_name = kwargs['det_name']
     my_d = bdv.Detector(det_name)
     
@@ -95,9 +98,3 @@ def job_main(kwargs):
     batch_loop(my_d, my_f, my_g4p, amplifier, g4_seed, total_events, instance_number)
     del my_g4p
 
-def main(kwargs):
-    scan_number = kwargs['scan']
-    for i in range(scan_number):
-        command = ' '.join(['python3', 'raser', '-b', 'gen_signal', '--job', str(i)] + sys.argv[3:]) # 'raser', '-sh', 'gen_signal'
-        print(command)
-        subprocess.run([command], shell=True)
\ No newline at end of file

particle/carrier_list.py

@@ -71,4 +71,101 @@ class CarrierListFromG4P:
         self.tracks_step = my_g4p.energy_steps[j]
         self.tracks_t_energy_deposition = my_g4p.edep_devices[j] #为什么不使用？
         self.ionized_pairs = [step*1e6/self.energy_loss for step in self.tracks_step]
-    
\ No newline at end of file
+
+class StripCarrierListFromG4P:
+    def __init__(self, material, my_g4p, batch):
+        if (material == "SiC"):
+            self.energy_loss = 8.4 #ev
+        elif (material == "Si"):
+            self.energy_loss = 3.6 #ev
+
+        if batch == 0:
+            h1 = ROOT.TH1F("Edep_device", "Energy deposition in Detector", 100, 0, max(my_g4p.edep_devices)*1.1)
+            for i in range (len(my_g4p.edep_devices)):
+                h1.Fill(my_g4p.edep_devices[i])
+            max_event_bin=h1.GetMaximumBin()
+            bin_wide=max(my_g4p.edep_devices)*1.1/100
+            for j in range (len(my_g4p.edep_devices)):
+                #compare to experimental data
+                if (my_g4p.edep_devices[j]<0.084 and my_g4p.edep_devices[j]>0.083):
+                    try_p=1
+                    for single_step in my_g4p.p_steps_current[j]:
+                        if abs(single_step[0]-my_g4p.p_steps_current[j][0][0])>5:
+                            try_p=0
+                    if try_p==1:
+                        self.batch_def(my_g4p,j)
+                        batch = 1
+                        break
+
+            if batch == 0:
+                print("the sensor didn't have particles hitted")
+                raise ValueError
+        else:
+            self.batch_def(my_g4p,batch)
+
+    def batch_def(self,my_g4p,j):
+        self.beam_number = j
+        self.track_position = [[single_step[0],single_step[1],single_step[2],1e-9] for single_step in my_g4p.p_steps_current[j]]
+        self.tracks_step = my_g4p.energy_steps[j]
+        self.tracks_t_energy_deposition = my_g4p.edep_devices[j] #为什么不使用？
+        self.ionized_pairs = [step*1e6/self.energy_loss for step in self.tracks_step]
+
+class PixelCarrierListFromG4P:
+    def __init__(self, my_d,my_g4p):
+        """
+        Description:
+            Events position and energy depositon
+        Parameters:
+            material : string
+                deciding the energy loss of MIP
+            my_g4p : Particles
+            batch : int
+                batch = 0: Single event, select particle with long enough track
+                batch != 0: Multi event, assign particle with batch number
+        Modify:
+            2022/10/25
+        """
+        batch = len(my_g4p.localposition)
+        layer = len(my_d.lt_z)
+        material = my_d.material
+        self.pixelsize_x = my_d.p_x
+        self.pixelsize_y = my_d.p_y
+        self.pixelsize_z = my_d.p_z
+        
+        if (material == "SiC"):
+            self.energy_loss = 8.4 #ev
+        elif (material == "Si"):
+            self.energy_loss = 3.6 #ev
+        
+        self.track_position, self.ionized_pairs= [],[]
+        self.layer= layer
+        for j in range(batch):
+            self.single_event(my_g4p,j)
+
+    def single_event(self,my_g4p,j):
+        s_track_position,s_energy= [],[]
+        for i in range(self.layer):
+            position = []
+            energy = []
+            name = "Layer_"+str(i)
+            #print(name)
+            for k in range(len(my_g4p.devicenames[j])):
+                px,py,pz = self.split_name(my_g4p.devicenames[j][k])
+                if name in my_g4p.devicenames[j][k]:
+                    tp = [0 for i in range(3)]
+                    tp[0] = my_g4p.localposition[j][k][0]+(px-0.5)*self.pixelsize_x
+                    tp[1] = my_g4p.localposition[j][k][1]+(py-0.5)*self.pixelsize_y
+                    tp[2] = my_g4p.localposition[j][k][2]+self.pixelsize_z/2
+                    position.append(tp) 
+                    energy.append(my_g4p.energy_steps[j][k])
+            s_track_position.append(position)
+            pairs = [step*1e6/self.energy_loss for step in energy]
+            s_energy.append(pairs)
+            del position,energy
+        self.track_position.append(s_track_position)
+        self.ionized_pairs.append(s_energy)
+        
+    def split_name(self,volume_name):
+        parts = volume_name.split('_')
+        return int(parts[1]),int(parts[2]),int(parts[4])
+

particle/g4_si_itk.py

@@ -511,7 +511,7 @@ def get_energy_deposition(my_g4p):
     create_path(path) 
     number = array('d',[999.])
     hittotal = array('d',[999.])
-    # number[0] = int(-ele_current.max_BB_height/18.8)
+    # number[0] = int(-ele_current.max_Broad_Band_height/18.8)
     # hittotal[0]=my_g4p.hittotal
     # fout = ROOT.TFile(path + "SiITk.root", "RECREATE")
     # t_out = ROOT.TTree("tree", "beam_number")

particle/g4_sic_energy_deposition.py

@@ -172,7 +172,7 @@ class SiC_LGAD_DetectorConstruction(g4b.G4VUserDetectorConstruction):#是g4b.G4V
         self.fMagFieldMessenger = g4b.G4GlobalMagFieldMessenger(fieldValue)
         self.fMagFieldMessenger.SetVerboseLevel(1)
 #创建粒子初始化信息，包括位置，动量，能量，
-class SIC_LGAD_PrimaryGeneratorAction(g4b.G4VUserPrimaryGeneratorAction):
+class SiC_LGAD_PrimaryGeneratorAction(g4b.G4VUserPrimaryGeneratorAction):
 
     def __init__(self):
         super().__init__()
@@ -294,7 +294,7 @@ class SiC_LGAD_aActionInitialization(g4b.G4VUserActionInitialization):
         self.SetUserAction(SiC_LGAD_RunAction())
 
     def Build(self):
-        self.SetUserAction(SIC_LGAD_PrimaryGeneratorAction())
+        self.SetUserAction(SiC_LGAD_PrimaryGeneratorAction())
         self.SetUserAction(SiC_LGAD_RunAction())
         eventAction = SiC_LGAD_aEventAction()
         self.SetUserAction(eventAction)