From 64c3f436dc5c25f84b7cb84110de703b205fa403 Mon Sep 17 00:00:00 2001
From: zhaosen <tokyosenzhao@163.com>
Date: Sat, 5 Oct 2024 19:16:48 +0800
Subject: [PATCH] half of gen_signal
---
current/cal_current.py | 2 +-
current/model.py | 2 +-
field/build_device.py | 1 +
field/devsim_field.py | 131 +++++++++++++++++-----------------
field/save_milestone.py | 15 ++--
field/solver_section.py | 12 ++--
gen_signal/draw_save.py | 3 +-
gen_signal/gen_signal_main.py | 14 ++--
particle/carrier_list.py | 12 ++--
9 files changed, 98 insertions(+), 94 deletions(-)
diff --git a/current/cal_current.py b/current/cal_current.py
index 6f947e6..c6948a2 100644
--- a/current/cal_current.py
+++ b/current/cal_current.py
@@ -20,7 +20,7 @@ from particle.carrier_list import CarrierListFromG4P, StripCarrierListFromG4P
from util.math import Vector, signal_convolution
from util.output import create_path
-t_bin = 50e-12
+t_bin = 100e-12
t_end = 10e-9
t_start = 0
delta_t = 10e-12
diff --git a/current/model.py b/current/model.py
index 5845768..49ac8f7 100644
--- a/current/model.py
+++ b/current/model.py
@@ -30,7 +30,7 @@ class Material:
self.hole_mass = 0.386*m_e
self.electron_mass = 0.26*m_e
self.avalanche_model = "vanOverstraeten"
- self.mobility_model = 'Reggiani'
+ self.mobility_model = 'Selberherr'
if self.mat_name == "SiC":
# 4H-SiC, use the longitudinal effective mass
self.permittivity=9.76
diff --git a/field/build_device.py b/field/build_device.py
index 56276a5..b684f01 100644
--- a/field/build_device.py
+++ b/field/build_device.py
@@ -50,6 +50,7 @@ class Detector:
self.absorber = self.device_dict['absorber']
self.amplifier = self.device_dict['amplifier']
+ self.read_out_contact = self.device_dict["read_out_contact"]
if "lgad3D" in self.det_model:
self.avalanche_bond = self.device_dict['avalanche_bond']
diff --git a/field/devsim_field.py b/field/devsim_field.py
index 679b6f1..f3f25d6 100644
--- a/field/devsim_field.py
+++ b/field/devsim_field.py
@@ -6,7 +6,12 @@
@Author : Henry Stone
@Version : 1.0
'''
-
+"""
+@File : Field_to_G4.py
+@Time : 2024/09/30
+@Author : Sen Zhao
+@Version : 2.0
+"""
import pickle
import ROOT
import numpy as np
@@ -16,23 +21,24 @@ from util.math import *
diff_res = 1e-5 # difference resolution in cm
class DevsimField:
- def __init__(self, device_name, dimension, voltage, read_ele_num, l_z):
+ def __init__(self, device_name, dimension, voltage, contact,read_ele_num):
self.name = device_name
self.voltage = voltage # float
self.dimension = dimension
self.read_ele_num = int(read_ele_num)
- self.l_z = l_z # used for planar weighting field TODO: auto weighting field
- DopingFile = "./output/field/{}/NetDoping_0.0V.pkl".format(self.name)
+
+
+ DopingFile = "./output/field/{}/NetDoping_0V.pkl".format(self.name)
PotentialFile = "./output/field/{}/Potential_{}V.pkl".format(self.name, self.voltage)
TrappingRate_pFile = "./output/field/{}/TrappingRate_p_{}V.pkl".format(self.name, self.voltage)
TrappingRate_nFile = "./output/field/{}/TrappingRate_n_{}V.pkl".format(self.name, self.voltage)
-
+ Weighting_Potential = "./output/field/{}/weightingfield/{}/Potential_{}V.pkl".format(self.name,contact, 1)
self.set_doping(DopingFile) #self.Doping
self.set_potential(PotentialFile) #self.Potential, self.x_efield, self.y_efield, self.z_efield
self.set_trap_p(TrappingRate_pFile) # self.TrappingRate_p
self.set_trap_n(TrappingRate_nFile) # self.TrappingRate_n
- self.set_w_p() #self.weighting_potential[]
+ self.set_w_p(Weighting_Potential) #self.weighting_potential[]
def set_doping(self, DopingFile):
try:
@@ -78,18 +84,30 @@ class DevsimField:
self.Potential = PotentialUniform
- def set_w_p(self):
- self.WeightingPotentials = [] #length = ele_num
- if self.read_ele_num == 1:
- print("Linear weighting potential loaded")
- pass
- elif self.read_ele_num >= 2:
- for i in range(self.read_ele_num):
- self.WeightingPotentials.append(strip_w_p(i))
- print("Weighting potential loaded for {}, strip {}".format(self.name, i+1))
- else:
- raise ValueError(self.read_ele_num)
+ def set_w_p(self,Weighting_PotentialFile):
+ try:
+ with open(Weighting_PotentialFile,'rb') as file:
+ Weighting_PotentialNotUniform=pickle.load(file)
+ print("Weighting_Potential file loaded for {}".format(self.name))
+ if Weighting_PotentialNotUniform['metadata']['dimension'] < self.dimension:
+ print("Weighting_Potential dimension not match")
+ return
+ except FileNotFoundError:
+ print("Weighting_Potential file not found, please run field simulation first")
+ print("or manually set the Weighting_Potential file")
+ return
+
+ if Weighting_PotentialNotUniform['metadata']['dimension'] == 1:
+ weighting_PotentialUniform = get_common_interpolate_1d(Weighting_PotentialNotUniform)
+ elif Weighting_PotentialNotUniform['metadata']['dimension'] == 2:
+ weighting_PotentialUniform = get_common_interpolate_2d(Weighting_PotentialNotUniform)
+ elif Weighting_PotentialNotUniform['metadata']['dimension'] == 3:
+ weighting_PotentialUniform = get_common_interpolate_3d(Weighting_PotentialNotUniform)
+
+ self.weighting_Potential = weighting_PotentialUniform
+
+
def set_trap_p(self, TrappingRate_pFile):
try:
with open(TrappingRate_pFile,'rb') as file:
@@ -134,6 +152,7 @@ class DevsimField:
self.TrappingRate_n = TrappingRate_nUniform
+
# DEVSIM dimension order: x, y, z
# RASER dimension order: z, x, y
@@ -203,7 +222,23 @@ class DevsimField:
E_x = - ((self.Potential(z, x) - self.Potential(z, x-diff_res))) / diff_res
except ValueError:
raise ValueError("Point {} might be out of bound x".format(x))
- return (E_x, 0, E_z)
+ try:
+ E_y = - ((self.Potential(z, x+diff_res/2) - self.Potential(z, x-diff_res/2))) / diff_res
+ except ValueError:
+ try:
+ E_y = - ((self.Potential(z, x+diff_res) - self.Potential(z, x))) / diff_res
+ except ValueError:
+ try:
+ E_y = - ((self.Potential(z, x) - self.Potential(z, x-diff_res))) / diff_res
+ except ValueError:
+ raise ValueError("Point {} might be out of bound x".format(y))
+ try:
+ return (E_x, 0, E_z)
+ except AttributeError:
+ try:
+ return (E_x,E_y,0)
+ except AttributeError:
+ return (0,E_y,E_z)
elif self.dimension == 3:
try:
@@ -238,15 +273,15 @@ class DevsimField:
raise ValueError("Point {} might be out of bound y".format(y))
return (E_x, E_y, E_z)
- def get_w_p(self, x, y, z, i):
- '''
- input: position in um
- output: weighting potential in 1
- '''
- if self.read_ele_num == 1:
- return linear_w_p(z, self.l_z)
- elif self.read_ele_num > 1:
- return self.WeightingPotentials[i].Interpolate(z, x)
+ def get_w_p(self, x, y, z,i): # used in cal current
+ x, y, z = x/1e4, y/1e4, z/1e4 # um to cm
+ if self.dimension == 1:
+ return self.weighting_Potential(z)
+ elif self.dimension == 2:
+ return self.weighting_Potential(z, x)
+ elif self.dimension == 3:
+ return self.weighting_Potential(z, x, y)
+
def get_trap_e(self, x, y, z):
'''
@@ -256,8 +291,10 @@ class DevsimField:
x, y, z = x/1e4, y/1e4, z/1e4 # um to cm
if self.dimension == 1:
return self.TrappingRate_n(z)
+
elif self.dimension == 2:
return self.TrappingRate_n(z, x)
+
elif self.dimension == 3:
return self.TrappingRate_n(z, x, y)
@@ -275,46 +312,6 @@ class DevsimField:
return self.TrappingRate_p(z, x, y)
-def linear_w_p(z, l_z):
- if z >= l_z:
- w_potential = 0
- elif z >= 1:
- w_potential = 1 - (1/(l_z-1)) * (z-1)
- else:
- w_potential = 1
- return w_potential
-
-def strip_w_p(ele_number):
- nx = 51
- ny = 321
- xmin, xmax = 0.0, 50.0
- ymin, ymax = 0.0, 320.0
- dx = (xmax - xmin) / (nx - 1)
- dy = (ymax - ymin) / (ny - 1)
-
- u = np.zeros((ny, nx))
- u[ele_number*75:(ele_number*75+20), 0] = 1.0
- u[:, -1] = 0.0
-
- max_iter = 100000
- tolerance = 1e-6
- for iteration in range(max_iter):
- u_old = u.copy()
- for i in range(1, ny - 1):
- for j in range(1, nx - 1):
- u[i, j] = (u[i+1, j] + u[i-1, j] + u[i, j+1] + u[i, j-1]) / 4
- diff = np.abs(u - u_old).max()
- if diff < tolerance:
- break
-
- x = np.linspace(xmin, xmax, nx)
- y = np.linspace(ymin, ymax, ny)
- w_potential=ROOT.TGraph2D()
- for i in range(len(y)):
- for j in range(len(x)):
- w_potential.SetPoint(int(i*len(x)+j),x[j]*6,y[i],u[i][j])
- return w_potential
-
if __name__ == "__main__":
testField = DevsimField("ITk-Si-strip", 2, -500.0, 4)
print(testField.get_e_field(100,100,50))
\ No newline at end of file
diff --git a/field/save_milestone.py b/field/save_milestone.py
index ed547d4..340372c 100644
--- a/field/save_milestone.py
+++ b/field/save_milestone.py
@@ -49,6 +49,7 @@ def milestone_save_2D(device, region, v, path):
Potential = np.array(devsim.get_node_model_values(device=device, region=region, name="Potential")) # get the potential data
TrappingRate_n = np.array(devsim.get_node_model_values(device=device, region=region, name="TrappingRate_n"))
TrappingRate_p = np.array(devsim.get_node_model_values(device=device, region=region, name="TrappingRate_p"))
+ NetDoping= np.array(devsim.get_node_model_values(device=device, region=region, name="NetDoping"))
devsim.element_from_edge_model(edge_model="ElectricField", device=device, region=region)
devsim.edge_average_model(device=device, region=region, node_model="x", edge_model="xmid")
@@ -70,8 +71,8 @@ def milestone_save_2D(device, region, v, path):
metadata['dimension'] = 2
names = ['Potential', 'TrappingRate_p', 'TrappingRate_n']
- # if v == 0:
- # names.append('NetDoping')
+ if v == 0:
+ names.append('NetDoping')
for name in names: # scalar field on mesh point (instead of on edge)
with open(os.path.join(path, "{}_{}V.pkl".format(name,v)),'wb') as file:
@@ -105,9 +106,9 @@ def milestone_save_wf_2D(device, region, v, path,contact):
metadata = {}
metadata['voltage'] = v
metadata['dimension'] = 2
-
+ save_wf_path = os.path.join(path,contact)
for name in ['Potential']: # scalar field on mesh point (instead of on edge)
- with open(os.path.join(path, "{}_{}_{}V.pkl".format(name,v,contact)),'wb') as file:
+ with open(os.path.join(save_wf_path, "{}_{}V.pkl".format(name,v)),'wb') as file:
data = {}
data['values'] = eval(name) # refer to the object with given name
merged_list = [x, y]
@@ -122,7 +123,7 @@ def milestone_save_3D(device, region, v, path):
Potential = np.array(devsim.get_node_model_values(device=device, region=region, name="Potential")) # get the potential data
TrappingRate_n = np.array(devsim.get_node_model_values(device=device, region=region, name="TrappingRate_n"))
TrappingRate_p = np.array(devsim.get_node_model_values(device=device, region=region, name="TrappingRate_p"))
-
+ NetDoping= np.array(devsim.get_node_model_values(device=device, region=region, name="NetDoping"))
devsim.element_from_edge_model(edge_model="ElectricField", device=device, region=region)
devsim.edge_average_model(device=device, region=region, node_model="x", edge_model="xmid")
devsim.edge_average_model(device=device, region=region, node_model="y", edge_model="ymid")
@@ -143,8 +144,8 @@ def milestone_save_3D(device, region, v, path):
metadata['dimension'] = 2
names = ['Potential', 'TrappingRate_p', 'TrappingRate_n']
- # if v == 0:
- # names.append('NetDoping')
+ if v == 0:
+ names.append('NetDoping')
for name in names: # scalar field on mesh point (instead of on edge)
with open(os.path.join(path, "{}_{}V.pkl".format(name,v)),'wb') as file:
diff --git a/field/solver_section.py b/field/solver_section.py
index 3bea26c..99d275b 100644
--- a/field/solver_section.py
+++ b/field/solver_section.py
@@ -195,7 +195,9 @@ def main (kwargs):
v_current=1
print("=======RASER info========\nBegin simulation WeightingField\n======================")
for contact in circuit_contacts:
- folder_path = os.path.join(path,"/weghtingfield/{}_{}_device".format(contact,v_current))
+ print(path)
+ folder_path = os.path.join(path, "weightingfield")
+ print(folder_path)
if not os.path.exists(folder_path):
os.makedirs(folder_path)
@@ -205,10 +207,10 @@ def main (kwargs):
value=0.0, acreal=paras['acreal'], acimag=paras['acimag'])
loop.initial_solver(contact=contact,set_contact_type=None,irradiation_label=irradiation_label,irradiation_flux=irradiation_flux,impact_label=impact_label)
- loop.loop_solver(circuit_contact=contact,v_current=v_current,area_factor=paras["area_factor"],path=path,device=device,region=region)
- if(paras['milestone_mode']==True and v_current%paras['milestone_step']==0.0):
- save_milestone.save_milestone(device=device, region=region, v=v_current, path=path,dimension=default_dimension,contact=contact,is_wf=is_wf)
- devsim.write_devices(file=os.path.join(folder_path,"weightingfield.dat"), type="tecplot")
+ loop.loop_solver(circuit_contact=contact,v_current=v_current,area_factor=paras["area_factor"],path=folder_path,device=device,region=region)
+
+ save_milestone.save_milestone(device=device, region=region, v=v_current, path=folder_path,dimension=default_dimension,contact=contact,is_wf=is_wf)
+ devsim.write_devices(file=os.path.join(folder_path,"weightingfield.dat"), type="tecplot")
elif is_wf == False:
v_current = 0
loop.initial_solver(contact=circuit_contacts,set_contact_type=None,irradiation_label=irradiation_label,irradiation_flux=irradiation_flux,impact_label=impact_label)
diff --git a/gen_signal/draw_save.py b/gen_signal/draw_save.py
index 06b38d5..e7e97a4 100644
--- a/gen_signal/draw_save.py
+++ b/gen_signal/draw_save.py
@@ -226,6 +226,7 @@ def draw_current(my_d, my_current, ele_current, read_ele_num, model, path, tag="
if ele_current[read_ele_num].GetMinimum() < 0:
rightmax = 1.1*ele_current[read_ele_num].GetMinimum()
+
else:
rightmax = 1.1*ele_current[read_ele_num].GetMaximum()
if rightmax == 0:
@@ -282,7 +283,7 @@ def cce(my_d,my_f,my_current, path):
for j in range(my_current.n_bin):
sum_charge=sum_charge+my_current.sum_cu[i].GetBinContent(j)*my_current.t_bin
charge.append(sum_charge/1.6e-19)
- print(charge)
+ print("===========RASER info================\nCollected Charge is {} C\n==============Result==============".format(list(charge)))
n=int(len(charge))
c1=ROOT.TCanvas("c1","canvas1",1000,1000)
cce=ROOT.TGraph(n,x,charge)
diff --git a/gen_signal/gen_signal_main.py b/gen_signal/gen_signal_main.py
index b015b18..aaab285 100644
--- a/gen_signal/gen_signal_main.py
+++ b/gen_signal/gen_signal_main.py
@@ -21,6 +21,7 @@ from elec import readout as rdo
from elec import ngspice_set_input as ngsip
from elec import ngspice as ng
from elec.set_pwl_input import set_pwl_input as pwlin
+import geant4_pybind as g4b
from .draw_save import energy_deposition, draw_drift_path, draw_current, cce
from util.output import output
@@ -50,12 +51,11 @@ def main(kwargs):
start = time.time()
det_name = kwargs['det_name']
- my_d = bdv.Detector(det_name)
-
+ my_d = bdv.Detector(det_name,devsim_solve_paras = None)
if kwargs['voltage'] != None:
- voltage = float(kwargs['voltage'])
+ voltage = str(kwargs['voltage'])
else:
- voltage = float(my_d.voltage)
+ voltage = str(my_d.voltage)
if kwargs['absorber'] != None:
absorber = kwargs['absorber']
@@ -67,12 +67,14 @@ def main(kwargs):
else:
amplifier = my_d.amplifier
- my_f = devfield.DevsimField(my_d.device, my_d.dimension, voltage, my_d.read_ele_num, my_d.l_z)
+ contact = my_d.read_out_contact
+
+ my_f = devfield.DevsimField(my_d.device, my_d.dimension, voltage,contact, my_d.read_ele_num)
g4_seed = random.randint(0,1e7)
my_g4p = g4t.Particles(my_d, absorber, g4_seed)
- if "strip" in det_name:
+ 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)
diff --git a/particle/carrier_list.py b/particle/carrier_list.py
index 783ccc3..e6e30a4 100644
--- a/particle/carrier_list.py
+++ b/particle/carrier_list.py
@@ -33,11 +33,11 @@ class CarrierListFromG4P:
if particle_number > 0:
batch=1
- if batch == 0:
- print("the sensor didn't have particles hitted")
- raise ValueError
- else:
- self.batch_def(my_g4p,batch)
+ if batch == 0:
+ print("=========RASER info ===========\nGeant4:the sensor didn't have particles hitted\n==========================")
+ raise ValueError
+ else:
+ self.batch_def(my_g4p,batch)
def batch_def(self,my_g4p,j):
self.beam_number = j
@@ -73,7 +73,7 @@ class StripCarrierListFromG4P:
break
if batch == 0:
- print("the sensor didn't have particles hitted")
+ print("=========RASER info ===========\nGeant4:the sensor didn't have particles hitted\n==========================")
raise ValueError
else:
self.batch_def(my_g4p,batch)
--
GitLab