diff --git a/__main__.py b/__main__.py
index 75b3acf88860e5c9032f5f78e3a578ff3ff48401..00b81dc1a736b7ccee38de73c6caee41c2dbaf50 100755
--- a/__main__.py
+++ b/__main__.py
@@ -43,6 +43,7 @@ parser_fpga.add_argument('label', help='LABEL to identify FPGA design')
parser_gen_signal = subparsers.add_parser('gen_signal', help='generate signal')
parser_gen_signal.add_argument('det_name', help='name of the detector')
+parser_gen_signal.add_argument('-l','--label', help='LABEL to identify signal generation method', default='signal')
parser_gen_signal.add_argument('-vol', '--voltage', type=str, help='bias voltage')
parser_gen_signal.add_argument('-abs', '--absorber', type=str, help='model of particle energy absorber')
parser_gen_signal.add_argument('-amp', '--amplifier', type=str, help='amplifier')
diff --git a/control/kei_2400b_control.py b/control/kei_2400b_bias_control.py
similarity index 100%
rename from control/kei_2400b_control.py
rename to control/kei_2400b_bias_control.py
diff --git a/control/kei_2400c_control.py b/control/kei_2400c_bias_control.py
similarity index 100%
rename from control/kei_2400c_control.py
rename to control/kei_2400c_bias_control.py
diff --git a/control/keithley2470ccontrol_hyj.py b/control/kei_2470c_bias_control.py
similarity index 100%
rename from control/keithley2470ccontrol_hyj.py
rename to control/kei_2470c_bias_control.py
diff --git a/control/keye4980acontrol.py b/control/key_e4980a_bias_control.py
similarity index 92%
rename from control/keye4980acontrol.py
rename to control/key_e4980a_bias_control.py
index ac040eab10b9522f026d6eed5fe3959d791b0c12..bf3dc11142b171608b97db210406d23efe661c5b 100644
--- a/control/keye4980acontrol.py
+++ b/control/key_e4980a_bias_control.py
@@ -2,7 +2,7 @@ import pyvisa as visa
import time
import warnings
-class keysighte4980a:
+class keysight_e4980a:
def __init__(self,resource_name):
instlist=visa.ResourceManager()
print(instlist.list_resources())
@@ -36,6 +36,6 @@ class keysighte4980a:
return float(cap)
if __name__=="__main__":
- lcr=keysighte4980a("USB0::0x2A8D::0x2F01::MY46516486::INSTR")
+ lcr=keysight_e4980a("USB0::0x2A8D::0x2F01::MY46516486::INSTR")
lcr.testIO()
lcr.set_trigger_internal()
diff --git a/control/mdo_3034_control.py b/control/mdo_3034_scope_control.py
similarity index 100%
rename from control/mdo_3034_control.py
rename to control/mdo_3034_scope_control.py
diff --git a/control/scan_amp.py b/control/scan_amp_2026.py
similarity index 91%
rename from control/scan_amp.py
rename to control/scan_amp_2026.py
index e8428c9107dd9c95917906096aff506c88ebc0bd..87b998f5965072259fd0ed8bafe2880c9bdb777d 100644
--- a/control/scan_amp.py
+++ b/control/scan_amp_2026.py
@@ -1,5 +1,5 @@
import raser.control.kei_2026_bias_control as kei2026
-import raser.control.tek_tds5054b_control as tds5054b
+import raser.control.tek_tds5054b_scope_control as tds5054b
import visa
import time
import pylab
diff --git a/control/scan_cv_2470.py b/control/scan_cv_2470.py
index 26f29a6e5e39dbd8c4afdbb7282af857002028a4..48e441425d646d65150d266e735b9f25d77718e4 100644
--- a/control/scan_cv_2470.py
+++ b/control/scan_cv_2470.py
@@ -1,5 +1,5 @@
-import keye4980acontrol
-import keithley2470ccontrol_hyj as kei2470
+import raser.control.key_e4980a_bias_control as keye4980a
+import raser.control.kei_2470c_bias_control as kei2470
import csv
import numpy as np
import platform
@@ -23,7 +23,7 @@ if platform.python_version().startswith('2'):
sys.exit()
-lcr=keye4980acontrol.keysighte4980a("USB0::0x2A8D::0x2F01::MY46516486::INSTR")
+lcr=keye4980a.keysight_e4980a("USB0::0x2A8D::0x2F01::MY46516486::INSTR")
lcr.set_frequency("1kHz")
lcr.set_voltage_level(0.1)
diff --git a/control/scan_iv.py b/control/scan_iv_2400.py
similarity index 90%
rename from control/scan_iv.py
rename to control/scan_iv_2400.py
index 36e15f5702bf21c1b6c20c92ffcac13d15d4ea5c..9fc524c7af07e2a1e260dbb310f85d376513341e 100644
--- a/control/scan_iv.py
+++ b/control/scan_iv_2400.py
@@ -1,4 +1,4 @@
-import raser.control.Kei2400c_control as kei2400
+import raser.control.kei_2400c_bias_control as kei2400
import visa
import time
import pylab
diff --git a/control/scan_strip_iv.py b/control/scan_iv_2400_strip.py
similarity index 90%
rename from control/scan_strip_iv.py
rename to control/scan_iv_2400_strip.py
index 7aec31aadb074ac391a60e0df397413b8937a4c4..9f7d1b6694983c88b38add774d17f760a9e349ad 100644
--- a/control/scan_strip_iv.py
+++ b/control/scan_iv_2400_strip.py
@@ -1,4 +1,4 @@
-import raser.control.kei_2400b_control as kei2400
+import raser.control.kei_2400b_bias_control as kei2400
import visa
import time
import pylab
diff --git a/control/scaniv.py b/control/scan_iv_2470.py
similarity index 93%
rename from control/scaniv.py
rename to control/scan_iv_2470.py
index 56412f9969f37d676049bcab47b9102aff91fab8..a60d3e87d2cae17acea7c95aee838cb2d51dc40f 100644
--- a/control/scaniv.py
+++ b/control/scan_iv_2470.py
@@ -1,5 +1,5 @@
#import Kei2400CControl as kei2400
-import keithley2470ccontrol_hyj as kei2470
+import raser.control.kei_2470c_bias_control as kei2470
# import visa
import pyvisa as visa
import time
diff --git a/control/tek_tds5054b_control.py b/control/tek_tds5054b_scope_control.py
similarity index 100%
rename from control/tek_tds5054b_control.py
rename to control/tek_tds5054b_scope_control.py
diff --git a/current/cal_current.py b/current/cal_current.py
index 14b69df59d61afc58d93ce3b9b685493b841257e..92ad06e0c895870e3fc35a786df33a48528c48fb 100644
--- a/current/cal_current.py
+++ b/current/cal_current.py
@@ -16,14 +16,13 @@ import numpy as np
import ROOT
from .model import Material
-from .model import Vector
+from util.math import Vector
from util.output import create_path
t_bin = 50e-12
t_end = 10e-9
t_start = 0
delta_t = 10e-12
-pixel = 25 #um
min_intensity = 1 # V/cm
class Carrier:
@@ -55,13 +54,9 @@ class Carrier:
self.z = z_init
self.t = t_init
self.t_end = t_end
- self.pixel = pixel
self.path = [[x_init, y_init, z_init, t_init]]
self.signal = [[] for j in range(int(read_ele_num))]
self.end_condition = 0
- self.diffuse_end_condition = 0
- self.row=0
- self.column=0
self.cal_mobility = Material(material).cal_mobility
self.charge = charge
@@ -79,9 +74,7 @@ class Carrier:
e_field = my_f.get_e_field(self.x,self.y,self.z)
intensity = Vector(e_field[0],e_field[1],e_field[2]).get_length()
mobility = Material(my_d.material)
- #mu = mobility.cal_mobility(my_d.temperature, my_d.doping_function(self.z+delta_z), self.charge, average_intensity)
- mu = mobility.cal_mobility(my_d.temperature, 1e12, self.charge, intensity)
- # TODO: rebuild the doping function or admit this as an approximation
+ mu = mobility.cal_mobility(my_d.temperature, my_f.get_doping(self.x, self.y, self.z), self.charge, intensity)
velocity_vector = [e_field[0]*mu, e_field[1]*mu, e_field[2]*mu] # cm/s
if(intensity > min_intensity):
@@ -163,80 +156,6 @@ class Carrier:
elif (self.t > t_end):
self.end_condition = "time out"
return self.end_condition
-
- def diffuse_single_step(self,my_d,my_f):
- delta_t=t_bin
- #e_field = my_f.get_e_field(self.x,self.y,self.z)
- intensity = 0
-
- kboltz=8.617385e-5 #eV/K
- mobility = Material(my_d.material)
- mu = mobility.cal_mobility(my_d.temperature, my_d.doping_function(self.z), self.charge, intensity)
- diffusion = (2.0*kboltz*mu*my_d.temperature*delta_t)**0.5
- #diffusion = 0.0
- dif_x=random.gauss(0.0,diffusion)*1e4
- dif_y=random.gauss(0.0,diffusion)*1e4
- dif_z=0
-
- if((self.x+dif_x)>=my_d.l_x):
- self.x = my_d.l_x
- elif((self.x+dif_x)<0):
- self.x = 0
- else:
- self.x = self.x+dif_x
- # y axis
- if((self.y+dif_y)>=my_d.l_y):
- self.y = my_d.l_y
- elif((self.y+dif_y)<0):
- self.y = 0
- else:
- self.y = self.y+dif_y
- # z axis
- if((self.z+dif_z)>=my_d.l_z):
- self.z = my_d.l_z
- elif((self.z+dif_z)<0):
- self.z = 0
- else:
- self.z = self.z+dif_z
- #time
- self.t = self.t+delta_t
- #record
- self.path.append([self.x,self.y,self.z,self.t])
-
- def diffuse_end(self,my_f):
- if (self.z<=0):
- # self.end_condition = "out of bound"
- self.diffuse_end_condition = "collect"
- return self.diffuse_end_condition
-
- def diffuse_not_in_sensor(self,my_d):
- if (self.x<=0) or (self.x>=my_d.l_x)\
- or (self.y<=0) or (self.y>=my_d.l_y)\
- or (self.z>=my_d.l_z):
- self.diffuse_end_condition = "out of bound"
- mod_x = self.x % self.pixel
- mod_y = self.y % self.pixel
- if ((mod_x> 7.5) & (mod_x<17.5)) & ((mod_y> 7.5) & (mod_y<17.5)) \
- & (self.t <= self.t_end):
- self.diffuse_end_condition = "collect"
- return self.diffuse_end_condition
-
- '''
- if (self.z<= 0) or (self.t >= self.t_end):
- self.diffuse_end_condition = "collect"
- #print("diffuse end")
- return self.diffuse_end_condition
- '''
-
- def pixel_position(self,my_f,my_d):
- if self.diffuse_end_condition == "collect":
- self.row = self.x // self.pixel
- self.column = self.y // self.pixel
- else:
- self.row = -1
- self.column = -1
- return self.row,self.column,abs(self.charge)
-
class CalCurrent:
@@ -557,108 +476,19 @@ class CalCurrentGain(CalCurrent):
self.negative_cu[j].Add(test_n)
test_n.Reset()
+
class CalCurrentG4P(CalCurrent):
def __init__(self, my_d, my_f, my_g4p, batch):
G4P_carrier_list = CarrierListFromG4P(my_d.material, my_g4p, batch)
super().__init__(my_d, my_f, G4P_carrier_list.ionized_pairs, G4P_carrier_list.track_position)
+
class CalCurrentStrip(CalCurrent):
def __init__(self, my_d, my_f, my_g4p, batch):
G4P_carrier_list = StripCarrierListFromG4P(my_d.material, my_g4p, batch)
super().__init__(my_d, my_f, G4P_carrier_list.ionized_pairs, G4P_carrier_list.track_position)
-class CalCurrentPixel:
- """Calculation of diffusion electrons in pixel detector"""
- def __init__(self, my_d, my_f, my_g4p):
- batch = len(my_g4p.localposition)
- layer = len(my_d.lt_z)
- G4P_carrier_list = PixelCarrierListFromG4P(my_d, my_g4p)
- self.collected_charge=[] #temp paras don't save as self.
- self.sum_signal = []
- self.event = []
- for k in range(batch):
- l_dict = {}
- signal_charge = []
- for j in range(layer):
- self.electrons = []
- self.charge,self.collected_charge = [],[]#same like before
- self.row,self.column=[],[]
- Hit = {'index':[],'charge':[]}
- #print(len(G4P_carrier_list.ionized_pairs[k][j]))
- print("%f pairs of carriers are generated from G4 in event_ %d layer %d" %(sum(G4P_carrier_list.ionized_pairs[k][j]),k,j))
- #print(G4P_carrier_list.track_position[k][j])
- for i in range(len(G4P_carrier_list.track_position[k][j])):
- electron = Carrier(G4P_carrier_list.track_position[k][j][i][0],\
- G4P_carrier_list.track_position[k][j][i][1],\
- G4P_carrier_list.track_position[k][j][i][2],\
- 0,\
- -1*G4P_carrier_list.ionized_pairs[k][j][i],\
- my_d.material,\
- 1)
- if not electron.not_in_sensor(my_d):
- self.electrons.append(electron)
- self.diffuse_loop(my_d,my_f)
-
- Xbins=int(my_d.l_x // electron.pixel)
- Ybins=int(my_d.l_y // electron.pixel)
- Xup=my_d.l_x // electron.pixel
- Yup=my_d.l_y // electron.pixel
- test_charge = ROOT.TH2F("charge", "charge",Xbins, 0, Xup, Ybins, 0, Yup)
- for i in range(len(self.row)):
- #test_charge.SetBinContent(int(self.row[i]),int(self.column[i]),self.charge[i])
- test_charge.Fill(self.row[i],self.column[i],self.charge[i])
-
- sum_fired = ROOT.TH2F("charge", "Pixel Detector charge",Xbins, 0, Xup, Ybins, 0, Yup)
- sum_fired.Add(test_charge)
-
- self.sum_charge = ROOT.TH2F("charge", "Pixel Detector charge",Xbins, 0, Xup, Ybins, 0, Yup)
- self.sum_charge.Add(test_charge)
-
- test_charge.Reset
- collected_charge=self.pixel_charge(my_d,Xbins,Ybins)
- signal_charge.append(collected_charge)
-
- Hit["index"],Hit["charge"] = self.pixel_fired(sum_fired,Xbins,Ybins)
- l_dict[j] = Hit
- print("%f electrons are collected in event_ %d,layer %d" %(sum(self.charge),k,j))
- self.sum_signal.append(signal_charge)
- self.event.append(l_dict)
- #print(signal_charge)
- del signal_charge
- #print(self.sum_signal)
- #print(self.event)
-
- def diffuse_loop(self, my_d, my_f):
- for electron in self.electrons:
- while not electron.diffuse_not_in_sensor(my_d):
- electron.diffuse_single_step(my_d, my_f)
- electron.diffuse_end(my_f)
- x,y,charge_quantity = electron.pixel_position(my_f,my_d)
- if (x != -1)&(y != -1):
- self.row.append(x)
- self.column.append(y)
- self.charge.append(charge_quantity)
-
- def pixel_charge(self,my_d,Xbins,Ybins):
- for x in range(Xbins):
- for y in range(Ybins):
- charge =self.sum_charge.GetBinContent(x,y)
- if (charge>0.2):
- self.collected_charge.append([x,y,charge])
- return self.collected_charge
-
- def pixel_fired(self,tot,Xbins,Ybins):
- Hit = {'index':[],'charge':[]}
- for x in range(Xbins):
- for y in range(Ybins):
- charge =tot.GetBinContent(x,y)
- if (charge>0.2):
- Hit["index"].append([x,y])
- Hit["charge"].append(charge)
- return Hit["index"],Hit["charge"]
-
-
class CalCurrentLaser(CalCurrent):
def __init__(self, my_d, my_f, my_l):
super().__init__(my_d, my_f, my_l.ionized_pairs, my_l.track_position)
@@ -750,65 +580,6 @@ class CarrierListFromG4P:
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])
-
class StripCarrierListFromG4P:
def __init__(self, material, my_g4p, batch):
diff --git a/current/cal_current_diffuse.py b/current/cal_current_diffuse.py
new file mode 100644
index 0000000000000000000000000000000000000000..4a6368c2ca31d97c6c422fc888d29f08669821c5
--- /dev/null
+++ b/current/cal_current_diffuse.py
@@ -0,0 +1,289 @@
+# -*- encoding: utf-8 -*-
+
+'''
+Description:
+ Simulate e-h pairs diffusing and calculate induced current in MAPS
+@Date : 2021/09/02 14:01:46
+@Author : Yuhang Tan, Haobo Wang
+@version : 2.0
+'''
+
+import random
+import math
+import os
+
+import numpy as np
+import ROOT
+
+from .model import Material
+from util.math import Vector
+
+t_bin = 50e-12
+t_end = 10e-9
+t_start = 0
+delta_t = 10e-12
+pixel = 25 #um
+min_intensity = 1 # V/cm
+
+class Carrier:
+ """
+ Description:
+ Definition of carriers and the record of their movement
+ Parameters:
+ x_init, y_init, z_init, t_init : float
+ initial space and time coordinates in um and s
+ charge : float
+ a set of drifting carriers, absolute value for number, sign for charge
+ Attributes:
+ x, y, z, t : float
+ space and time coordinates in um and s
+ path : float[]
+ recording the carrier path in [x, y, z, t]
+ charge : float
+ a set of drifting carriers, absolute value for number, sign for charge
+ signal : float[]
+ the generated signal current on the reading electrode
+ end_condition : 0/string
+ tag of how the carrier ended drifting
+ Modify:
+ 2022/10/28
+ """
+ def __init__(self, x_init, y_init, z_init, t_init, charge, read_ele_num):
+ self.x = x_init
+ self.y = y_init
+ self.z = z_init
+ self.t = t_init
+ self.t_end = t_end
+ self.pixel = pixel
+ self.path = [[x_init, y_init, z_init, t_init]]
+ self.signal = [[] for j in range(int(read_ele_num))]
+ self.diffuse_end_condition = 0
+ self.row=0
+ self.column=0
+
+ self.charge = charge
+
+ def diffuse_single_step(self,my_d,my_f):
+ delta_t=t_bin
+ #e_field = my_f.get_e_field(self.x,self.y,self.z)
+ intensity = 0
+
+ kboltz=8.617385e-5 #eV/K
+ mobility = Material(my_d.material)
+ mu = mobility.cal_mobility(my_d.temperature, my_f.get_doping(self.x, self.y, self.z), self.charge, intensity)
+ diffusion = (2.0*kboltz*mu*my_d.temperature*delta_t)**0.5
+ #diffusion = 0.0
+ dif_x=random.gauss(0.0,diffusion)*1e4
+ dif_y=random.gauss(0.0,diffusion)*1e4
+ dif_z=0
+
+ if((self.x+dif_x)>=my_d.l_x):
+ self.x = my_d.l_x
+ elif((self.x+dif_x)<0):
+ self.x = 0
+ else:
+ self.x = self.x+dif_x
+ # y axis
+ if((self.y+dif_y)>=my_d.l_y):
+ self.y = my_d.l_y
+ elif((self.y+dif_y)<0):
+ self.y = 0
+ else:
+ self.y = self.y+dif_y
+ # z axis
+ if((self.z+dif_z)>=my_d.l_z):
+ self.z = my_d.l_z
+ elif((self.z+dif_z)<0):
+ self.z = 0
+ else:
+ self.z = self.z+dif_z
+ #time
+ self.t = self.t+delta_t
+ #record
+ self.path.append([self.x,self.y,self.z,self.t])
+
+ def diffuse_end(self,my_f):
+ if (self.z<=0):
+ # self.end_condition = "out of bound"
+ self.diffuse_end_condition = "collect"
+ return self.diffuse_end_condition
+
+ def diffuse_not_in_sensor(self,my_d):
+ if (self.x<=0) or (self.x>=my_d.l_x)\
+ or (self.y<=0) or (self.y>=my_d.l_y)\
+ or (self.z>=my_d.l_z):
+ self.diffuse_end_condition = "out of bound"
+ mod_x = self.x % self.pixel
+ mod_y = self.y % self.pixel
+ if ((mod_x> 7.5) & (mod_x<17.5)) & ((mod_y> 7.5) & (mod_y<17.5)) \
+ & (self.t <= self.t_end):
+ self.diffuse_end_condition = "collect"
+ return self.diffuse_end_condition
+
+ '''
+ if (self.z<= 0) or (self.t >= self.t_end):
+ self.diffuse_end_condition = "collect"
+ #print("diffuse end")
+ return self.diffuse_end_condition
+ '''
+
+ def pixel_position(self,my_f,my_d):
+ if self.diffuse_end_condition == "collect":
+ self.row = self.x // self.pixel
+ self.column = self.y // self.pixel
+ else:
+ self.row = -1
+ self.column = -1
+ return self.row,self.column,abs(self.charge)
+
+
+
+class CalCurrentPixel:
+ """Calculation of diffusion electrons in pixel detector"""
+ def __init__(self, my_d, my_f, my_g4p):
+ batch = len(my_g4p.localposition)
+ layer = len(my_d.lt_z)
+ G4P_carrier_list = PixelCarrierListFromG4P(my_d, my_g4p)
+ self.collected_charge=[] #temp paras don't save as self.
+ self.sum_signal = []
+ self.event = []
+ for k in range(batch):
+ l_dict = {}
+ signal_charge = []
+ for j in range(layer):
+ self.electrons = []
+ self.charge,self.collected_charge = [],[]#same like before
+ self.row,self.column=[],[]
+ Hit = {'index':[],'charge':[]}
+ #print(len(G4P_carrier_list.ionized_pairs[k][j]))
+ print("%f pairs of carriers are generated from G4 in event_ %d layer %d" %(sum(G4P_carrier_list.ionized_pairs[k][j]),k,j))
+ #print(G4P_carrier_list.track_position[k][j])
+ for i in range(len(G4P_carrier_list.track_position[k][j])):
+ electron = Carrier(G4P_carrier_list.track_position[k][j][i][0],\
+ G4P_carrier_list.track_position[k][j][i][1],\
+ G4P_carrier_list.track_position[k][j][i][2],\
+ 0,\
+ -1*G4P_carrier_list.ionized_pairs[k][j][i],\
+ 1)
+ if not electron.diffuse_not_in_sensor(my_d):
+ self.electrons.append(electron)
+ self.diffuse_loop(my_d,my_f)
+
+ Xbins=int(my_d.l_x // electron.pixel)
+ Ybins=int(my_d.l_y // electron.pixel)
+ Xup=my_d.l_x // electron.pixel
+ Yup=my_d.l_y // electron.pixel
+ test_charge = ROOT.TH2F("charge", "charge",Xbins, 0, Xup, Ybins, 0, Yup)
+ for i in range(len(self.row)):
+ #test_charge.SetBinContent(int(self.row[i]),int(self.column[i]),self.charge[i])
+ test_charge.Fill(self.row[i],self.column[i],self.charge[i])
+
+ sum_fired = ROOT.TH2F("charge", "Pixel Detector charge",Xbins, 0, Xup, Ybins, 0, Yup)
+ sum_fired.Add(test_charge)
+
+ self.sum_charge = ROOT.TH2F("charge", "Pixel Detector charge",Xbins, 0, Xup, Ybins, 0, Yup)
+ self.sum_charge.Add(test_charge)
+
+ test_charge.Reset
+ collected_charge=self.pixel_charge(my_d,Xbins,Ybins)
+ signal_charge.append(collected_charge)
+
+ Hit["index"],Hit["charge"] = self.pixel_fired(sum_fired,Xbins,Ybins)
+ l_dict[j] = Hit
+ print("%f electrons are collected in event_ %d,layer %d" %(sum(self.charge),k,j))
+ self.sum_signal.append(signal_charge)
+ self.event.append(l_dict)
+ #print(signal_charge)
+ del signal_charge
+ #print(self.sum_signal)
+ #print(self.event)
+
+ def diffuse_loop(self, my_d, my_f):
+ for electron in self.electrons:
+ while not electron.diffuse_not_in_sensor(my_d):
+ electron.diffuse_single_step(my_d, my_f)
+ electron.diffuse_end(my_f)
+ x,y,charge_quantity = electron.pixel_position(my_f,my_d)
+ if (x != -1)&(y != -1):
+ self.row.append(x)
+ self.column.append(y)
+ self.charge.append(charge_quantity)
+
+ def pixel_charge(self,my_d,Xbins,Ybins):
+ for x in range(Xbins):
+ for y in range(Ybins):
+ charge =self.sum_charge.GetBinContent(x,y)
+ if (charge>0.2):
+ self.collected_charge.append([x,y,charge])
+ return self.collected_charge
+
+ def pixel_fired(self,tot,Xbins,Ybins):
+ Hit = {'index':[],'charge':[]}
+ for x in range(Xbins):
+ for y in range(Ybins):
+ charge =tot.GetBinContent(x,y)
+ if (charge>0.2):
+ Hit["index"].append([x,y])
+ 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])
+
diff --git a/current/model.py b/current/model.py
index d14752479af598f3ba22813012a9da0dcf5412cc..339932b0084a7c4c1bf10d75a45e0f3c9eb7b8bd 100644
--- a/current/model.py
+++ b/current/model.py
@@ -330,38 +330,6 @@ class Material:
return coefficient
-class Vector:
- def __init__(self,a1,a2,a3):
- self.components = [a1,a2,a3]
-
- def cross(self,Vector_b):
- """ Get vector cross product of self and another Vector"""
- o1 = self.components[1]*Vector_b.components[2]-self.components[2]*Vector_b.components[1]
- o2 = self.components[2]*Vector_b.components[0]-self.components[0]*Vector_b.components[2]
- o3 = self.components[0]*Vector_b.components[1]-self.components[1]*Vector_b.components[0]
- return Vector(o1,o2,o3)
-
- def get_length(self):
- " Return length of self"
- return math.sqrt(self.components[0]*self.components[0]+self.components[1]*self.components[1]+self.components[2]*self.components[2])
-
- def add(self,Vector_b):
- " Return the sum of two Vectors. eg: [1,2,3]+[1,2,3] = [2,4,6]"
- o1 = self.components[0]+Vector_b.components[0]
- o2 = self.components[1]+Vector_b.components[1]
- o3 = self.components[2]+Vector_b.components[2]
- return Vector(o1,o2,o3)
-
- def sub(self,Vector_b):
- " Return the subtraction of two Vectors. eg: [1,2,3]-[1,2,3] = [0,0,0]"
- o1 = self.components[0]-Vector_b.components[0]
- o2 = self.components[1]-Vector_b.components[1]
- o3 = self.components[2]-Vector_b.components[2]
- return Vector(o1,o2,o3)
-
- def mul(self,k):
- " Return Vector multiplied by number. eg: 2 * [1,2,3] = [2,4,6]"
- return Vector(self.components[0]*k,self.components[1]*k,self.components[2]*k)
def main():
if not (os.path.exists("./output/model")):
diff --git a/draw/__init__.py b/draw/__init__.py
index 56d912c04af8f1cd7e0be6f7dfd6ecfaeaef12c4..59dc25f741ef22e6c47f5c9451ce8176a6edb393 100644
--- a/draw/__init__.py
+++ b/draw/__init__.py
@@ -5,15 +5,17 @@ from . import cv
from . import compare_iv
from . import compare_cv
from . import draw_iv_cv_paper9
-from . import alpha_ana_SICAR
+from . import alpha_ana_sicar
from . import data_convolution
+from . import draw_figure_sicar as dfc
+
def main(kwargs):
label = kwargs['label']
if label == 'sicar1.1.8':
input_dir = '/publicfs/atlas/atlasnew/silicondet/itk/raser/wangkeqi/sicar1.1.8'
output_dir = '/afs/ihep.ac.cn/users/w/wangkeqi/raser/output/fig'
- draw_figure(input_dir, output_dir, label)
+ dfc.draw(input_dir, output_dir, label)
elif label == 'compare_sicar1.1.8_iv':
iv.main(label)
elif label == 'compare_sicar1.1.8_cv':
@@ -21,22 +23,22 @@ def main(kwargs):
elif label == 'sicar1.1.8-1':
input_dir = '/publicfs/atlas/atlasnew/silicondet/itk/raser/wangkeqi/sicar1.1.8'
output_dir = '/afs/ihep.ac.cn/users/w/wangkeqi/raser/output/fig'
- draw_figure(input_dir, output_dir, label)
+ dfc.draw(input_dir, output_dir, label)
elif label == 'sicar1.1.8-2':
input_dir = '/publicfs/atlas/atlasnew/silicondet/itk/raser/wangkeqi/sicar1.1.8'
output_dir = '/afs/ihep.ac.cn/users/w/wangkeqi/raser/output/fig'
- draw_figure(input_dir, output_dir, label)
+ dfc.draw(input_dir, output_dir, label)
elif label == 'itk_md8_data_v1':
input_dir = '/publicfs/atlas/atlasnew/silicondet/itk/raser/lizhan/itkmd8/itkmd8data'
output_dir = '/afs/ihep.ac.cn/users/l/lizhan/disk/scrathfs/raser/output/fig'
- draw_figure(input_dir, output_dir, label,xtitle_iv="Reverse Bias Voltage [V]",ytitle_iv="Current [nA]",
+ dfc.draw(input_dir, output_dir, label,xtitle_iv="Reverse Bias Voltage [V]",ytitle_iv="Current [nA]",
xtitle_cv="Reverse Bias Voltage [V]",ytitle_cv="Capacitance [pF]",
xlowerlimit_iv=0,xupperlimit_iv=700,ylowerlimit_iv=1e-11,yupperlimit_iv=1e-5,ylogscale_iv=0,
xlowerlimit_cv=0,xupperlimit_cv=400,ylowerlimit_cv=0,yupperlimit_cv=1e2,ylogscale_cv=0)
elif label == 'itk_md8_sim_v1':
input_dir = '/publicfs/atlas/atlasnew/silicondet/itk/raser/lizhan/itkmd8/itkmd8sim'
output_dir = '/afs/ihep.ac.cn/users/l/lizhan/disk/scrathfs/raser/output/fig'
- draw_figure(input_dir, output_dir, label,xtitle_iv="Reverse Bias Voltage [V]",ytitle_iv="Current [nA]",
+ dfc.draw(input_dir, output_dir, label,xtitle_iv="Reverse Bias Voltage [V]",ytitle_iv="Current [nA]",
xtitle_cv="Reverse Bias Voltage [V]",ytitle_cv="Capacitance [pF]",
xlowerlimit_iv=0,xupperlimit_iv=700,ylowerlimit_iv=1e-11,yupperlimit_iv=1e-5,ylogscale_iv=0,
xlowerlimit_cv=0,xupperlimit_cv=400,ylowerlimit_cv=0,yupperlimit_cv=1e2,ylogscale_cv=0)
@@ -48,14 +50,14 @@ def main(kwargs):
elif label == 'itk_atlas18_sim_v1':
input_dir = '/publicfs/atlas/atlasnew/silicondet/itk/raser/lizhan/atlas18/sim'
output_dir = '/afs/ihep.ac.cn/users/l/lizhan/disk/scrathfs/raser/output/fig'
- draw_figure(input_dir, output_dir, label,xtitle_iv="Reverse Bias Voltage [V]",ytitle_iv="Current [A]",
+ dfc.draw(input_dir, output_dir, label,xtitle_iv="Reverse Bias Voltage [V]",ytitle_iv="Current [A]",
xtitle_cv="Reverse Bias Voltage [V]",ytitle_cv="Capacitance [pF]",
xlowerlimit_iv=0,xupperlimit_iv=700,ylowerlimit_iv=1e-11,yupperlimit_iv=1e-5,ylogscale_iv=0,
xlowerlimit_cv=0,xupperlimit_cv=400,ylowerlimit_cv=0,yupperlimit_cv=1e2,ylogscale_cv=0)
elif label == 'itk_atlas18_data_v1':
input_dir = '/publicfs/atlas/atlasnew/silicondet/itk/raser/lizhan/atlas18/data'
output_dir = '/afs/ihep.ac.cn/users/l/lizhan/disk/scrathfs/raser/output/fig'
- draw_figure(input_dir, output_dir, label,xtitle_iv="Reverse Bias Voltage [V]",ytitle_iv="Current [nA]",
+ dfc.draw(input_dir, output_dir, label,xtitle_iv="Reverse Bias Voltage [V]",ytitle_iv="Current [nA]",
xtitle_cv="Reverse Bias Voltage [V]",ytitle_cv="Capacitance [pF]",
xlowerlimit_iv=0,xupperlimit_iv=700,ylowerlimit_iv=1e-11,yupperlimit_iv=1e-5,ylogscale_iv=0,
xlowerlimit_cv=0,xupperlimit_cv=400,ylowerlimit_cv=0,yupperlimit_cv=1e2,ylogscale_cv=0)
@@ -92,117 +94,9 @@ def main(kwargs):
root_fullpath = os.path.join(csv_folder, root_filename)
draw_iv_cv_paper9.create_root_file(csv_filename, root_fullpath)
draw_iv_cv_paper9.main()
- alpha_ana_SICAR.main()
+ alpha_ana_sicar.main()
data_convolution.landau_mirror()
data_convolution.energy_sim()
else:
raise NameError(label)
-
-def draw_figure(input_dir, output_dir, label,xtitle_iv="Reverse Bias Voltage [V]",ytitle_iv="Current [A]",
- xtitle_cv="Reverse Bias Voltage [V]",ytitle_cv="Capacitance [pF]",
- xlowerlimit_iv=0,xupperlimit_iv=510,ylowerlimit_iv=1e-11,yupperlimit_iv=1e-5,ylogscale_iv=0,
- xlowerlimit_cv=0,xupperlimit_cv=399.99,ylowerlimit_cv=0,yupperlimit_cv=1e2,ylogscale_cv=0):
-
- com_name = []
- for file in os.listdir(input_dir):
- if file.endswith('.root'):
- com_name.append(file)
- for name in com_name:
- if label == 'sicar1.1.8' and not name.startswith('sicar1.1.8'):
- continue
- elif label == 'sicar1.1.8-1' and not name.startswith('sicar1.1.8-1_'):
- continue
- elif label == 'sicar1.1.8-2' and not name.startswith('sicar1.1.8-2_'):
- continue
- name = name.split('.root')[0]
-
- input_file = os.path.join(input_dir, name + '.root')
- output_file = os.path.join(output_dir, name + '.root')
- pdf_file = os.path.join(output_dir, name + '.pdf')
- png_file = os.path.join(output_dir, name + '.png')
-
- if name.endswith('iv'):
- file = ROOT.TFile(input_file, "READ")
- tree = file.Get("myTree")
- graph = ROOT.TGraph()
-
- for i, event in enumerate(tree):
- if label in ['itk_md8_data','itk_atlas18_data']:
- x = event.Voltage_V
- x = abs(x)
- y = event.Current_nA
- y = abs(y)*1e-9
- elif label == 'itk_atlas18_sim':
- x = event.Voltage
- x = abs(x)
- y = event.Current
- y = abs(y)*1e-9
- else:
- x = event.Value
- x = abs(x)
- y = event.Reading
- y = abs(y)
- graph.SetPoint(i, x, y)
-
- draw_with_options(graph,name,output_file,pdf_file,png_file,xtitle_iv,ytitle_iv,
- xlowerlimit_iv,xupperlimit_iv,ylowerlimit_iv,yupperlimit_iv,ylogscale_iv)
- #problem: unable to change y limits
-
- if name.endswith('cv'):
- file = ROOT.TFile(input_file, "READ")
- tree = file.Get("myTree")
- graph = ROOT.TGraph()
- for i, event in enumerate(tree):
- x = event.Voltage
- x = abs(x)
- y = event.Capacitance
- y = abs(y)
- graph.SetPoint(i, x, y)
-
- draw_with_options(graph,name,output_file,pdf_file,png_file,xtitle_cv,ytitle_cv,
- xlowerlimit_cv,xupperlimit_cv,ylowerlimit_cv,yupperlimit_cv,ylogscale_cv)
-
-
-def draw_with_options(graph,name,output_file,pdf_file,png_file,xtitle,ytitle,
- xlowerlimit,xupperlimit,ylowerlimit,yupperlimit,ylogscale):
- graph.SetNameTitle("")
- graph.SetLineWidth(1)
- graph.SetMarkerColor(ROOT.kBlack)
- graph.SetMarkerStyle(24)
- graph.SetMarkerSize(1)
-
- graph.GetXaxis().SetTitle(xtitle)
- graph.GetXaxis().SetLimits(xlowerlimit,xupperlimit)
- graph.GetXaxis().CenterTitle()
- graph.GetXaxis().SetTitleOffset(1.4)
- graph.GetXaxis().SetTitleSize(0.05)
- graph.GetXaxis().SetLabelSize(0.05)
- graph.GetXaxis().SetNdivisions(505)
-
- graph.GetYaxis().SetLimits(ylowerlimit,yupperlimit)
- graph.GetYaxis().SetTitle(ytitle)
- graph.GetYaxis().CenterTitle()
- graph.GetYaxis().SetTitleOffset(1.8)
- graph.GetYaxis().SetTitleSize(0.05)
- graph.GetYaxis().SetLabelSize(0.05)
- graph.Draw("AP")
-
- c = ROOT.TCanvas("c","c",500,500)
- c.SetLeftMargin(0.22)
- c.SetBottomMargin(0.16)
- legend = ROOT.TLegend(0.27,0.67,0.62,0.80)
- c.SetGrid()
- c.SetFrameLineWidth(5)
-
- legend.SetTextSize(0.04)
- legend.AddEntry(graph,name.split('_')[0])
-
- c.cd()
- c.SetLogy(ylogscale)
- graph.Draw()
- legend.Draw()
-
- c.SaveAs(output_file)
- c.SaveAs(pdf_file)
- c.SaveAs(png_file)
- del c
\ No newline at end of file
+
\ No newline at end of file
diff --git a/draw/alpha_ana_SICAR.py b/draw/alpha_ana_sicar.py
similarity index 100%
rename from draw/alpha_ana_SICAR.py
rename to draw/alpha_ana_sicar.py
diff --git a/draw/draw_figure_sicar.py b/draw/draw_figure_sicar.py
new file mode 100644
index 0000000000000000000000000000000000000000..ab74bcf6748272e4945e9152a00ee8702b10f8c5
--- /dev/null
+++ b/draw/draw_figure_sicar.py
@@ -0,0 +1,110 @@
+import os
+import ROOT
+def draw(input_dir, output_dir, label,xtitle_iv="Reverse Bias Voltage [V]",ytitle_iv="Current [A]",
+ xtitle_cv="Reverse Bias Voltage [V]",ytitle_cv="Capacitance [pF]",
+ xlowerlimit_iv=0,xupperlimit_iv=510,ylowerlimit_iv=1e-11,yupperlimit_iv=1e-5,ylogscale_iv=0,
+ xlowerlimit_cv=0,xupperlimit_cv=399.99,ylowerlimit_cv=0,yupperlimit_cv=1e2,ylogscale_cv=0):
+
+ com_name = []
+ for file in os.listdir(input_dir):
+ if file.endswith('.root'):
+ com_name.append(file)
+ for name in com_name:
+ if label == 'sicar1.1.8' and not name.startswith('sicar1.1.8'):
+ continue
+ elif label == 'sicar1.1.8-1' and not name.startswith('sicar1.1.8-1_'):
+ continue
+ elif label == 'sicar1.1.8-2' and not name.startswith('sicar1.1.8-2_'):
+ continue
+ name = name.split('.root')[0]
+
+ input_file = os.path.join(input_dir, name + '.root')
+ output_file = os.path.join(output_dir, name + '.root')
+ pdf_file = os.path.join(output_dir, name + '.pdf')
+ png_file = os.path.join(output_dir, name + '.png')
+
+ if name.endswith('iv'):
+ file = ROOT.TFile(input_file, "READ")
+ tree = file.Get("myTree")
+ graph = ROOT.TGraph()
+
+ for i, event in enumerate(tree):
+ if label in ['itk_md8_data','itk_atlas18_data']:
+ x = event.Voltage_V
+ x = abs(x)
+ y = event.Current_nA
+ y = abs(y)*1e-9
+ elif label == 'itk_atlas18_sim':
+ x = event.Voltage
+ x = abs(x)
+ y = event.Current
+ y = abs(y)*1e-9
+ else:
+ x = event.Value
+ x = abs(x)
+ y = event.Reading
+ y = abs(y)
+ graph.SetPoint(i, x, y)
+
+ draw_with_options(graph,name,output_file,pdf_file,png_file,xtitle_iv,ytitle_iv,
+ xlowerlimit_iv,xupperlimit_iv,ylowerlimit_iv,yupperlimit_iv,ylogscale_iv)
+ #problem: unable to change y limits
+
+ if name.endswith('cv'):
+ file = ROOT.TFile(input_file, "READ")
+ tree = file.Get("myTree")
+ graph = ROOT.TGraph()
+ for i, event in enumerate(tree):
+ x = event.Voltage
+ x = abs(x)
+ y = event.Capacitance
+ y = abs(y)
+ graph.SetPoint(i, x, y)
+
+ draw_with_options(graph,name,output_file,pdf_file,png_file,xtitle_cv,ytitle_cv,
+ xlowerlimit_cv,xupperlimit_cv,ylowerlimit_cv,yupperlimit_cv,ylogscale_cv)
+
+
+def draw_with_options(graph,name,output_file,pdf_file,png_file,xtitle,ytitle,
+ xlowerlimit,xupperlimit,ylowerlimit,yupperlimit,ylogscale):
+ graph.SetNameTitle("")
+ graph.SetLineWidth(1)
+ graph.SetMarkerColor(ROOT.kBlack)
+ graph.SetMarkerStyle(24)
+ graph.SetMarkerSize(1)
+
+ graph.GetXaxis().SetTitle(xtitle)
+ graph.GetXaxis().SetLimits(xlowerlimit,xupperlimit)
+ graph.GetXaxis().CenterTitle()
+ graph.GetXaxis().SetTitleOffset(1.4)
+ graph.GetXaxis().SetTitleSize(0.05)
+ graph.GetXaxis().SetLabelSize(0.05)
+ graph.GetXaxis().SetNdivisions(505)
+
+ graph.GetYaxis().SetLimits(ylowerlimit,yupperlimit)
+ graph.GetYaxis().SetTitle(ytitle)
+ graph.GetYaxis().CenterTitle()
+ graph.GetYaxis().SetTitleOffset(1.8)
+ graph.GetYaxis().SetTitleSize(0.05)
+ graph.GetYaxis().SetLabelSize(0.05)
+ graph.Draw("AP")
+
+ c = ROOT.TCanvas("c","c",500,500)
+ c.SetLeftMargin(0.22)
+ c.SetBottomMargin(0.16)
+ legend = ROOT.TLegend(0.27,0.67,0.62,0.80)
+ c.SetGrid()
+ c.SetFrameLineWidth(5)
+
+ legend.SetTextSize(0.04)
+ legend.AddEntry(graph,name.split('_')[0])
+
+ c.cd()
+ c.SetLogy(ylogscale)
+ graph.Draw()
+ legend.Draw()
+
+ c.SaveAs(output_file)
+ c.SaveAs(pdf_file)
+ c.SaveAs(png_file)
+ del c
\ No newline at end of file
diff --git a/field/build_device.py b/field/build_device.py
index 921d5fa4c7254f224bd388b405d676297605e74c..9bb4dffa6194b95626aea4caeae47aacf4582de5 100644
--- a/field/build_device.py
+++ b/field/build_device.py
@@ -1,14 +1,16 @@
#!/usr/bin/env python3
# -*- encoding: utf-8 -*-
-import devsim
-from . import model_create
-
-from util.output import output
import json
import os
+import devsim
import matplotlib.pyplot
+import numpy as np
+
+from . import model_create
+from util.output import output
+from util.math import *
class Detector:
"""
@@ -65,8 +67,10 @@ class Detector:
if "strip" in self.det_name or "Strip" in self.det_name:
# TODO: change this into model
self.read_ele_num = self.device_dict['read_ele_num']
+ else:
+ self.read_ele_num = 1
- if "pixeldetector" in self.det_model:
+ if "pixel" in self.det_model:
self.p_x = self.device_dict['px']
self.p_y = self.device_dict['py']
self.p_z = self.device_dict['pz']
@@ -148,13 +152,13 @@ class Detector:
elif self.control_dict["ac-weightfield"] == False:
pass
if 'Acceptors_ir' in self.device_dict['doping']:
- model_create.CreateNodeModel(self.device, self.region, "Acceptors", self.device_dict['doping']['Acceptors']+"+"+self.device_dict['doping']['Acceptors_ir'])
+ model_create.CreateNodeModel(self.device, self.region, "Acceptors", self.device_dict['doping']['Acceptors']+"+"+self.device_dict['doping']['Acceptors_ir'])
else:
- model_create.CreateNodeModel(self.device, self.region, "Acceptors", self.device_dict['doping']['Acceptors'])
+ model_create.CreateNodeModel(self.device, self.region, "Acceptors", self.device_dict['doping']['Acceptors'])
if 'Donors_ir' in self.device_dict['doping']:
- model_create.CreateNodeModel(self.device, self.region, "Donors", self.device_dict['doping']['Donors']+"+"+self.device_dict['doping']['Donors_ir'])
+ model_create.CreateNodeModel(self.device, self.region, "Donors", self.device_dict['doping']['Donors']+"+"+self.device_dict['doping']['Donors_ir'])
else:
- model_create.CreateNodeModel(self.device, self.region, "Donors", self.device_dict['doping']['Donors'])
+ model_create.CreateNodeModel(self.device, self.region, "Donors", self.device_dict['doping']['Donors'])
model_create.CreateNodeModel(self.device, self.region, "NetDoping", "Donors-Acceptors")
devsim.edge_from_node_model(device=self.device, region=self.region, node_model="Acceptors")
devsim.edge_from_node_model(device=self.device, region=self.region, node_model="NetDoping")
diff --git a/field/devsim_field.py b/field/devsim_field.py
index 80daf40f36709ec7950dc6a6b7f7f5612cc1e189..679b6f12a813e236009dc670bbbeea1e0d4e0b9c 100644
--- a/field/devsim_field.py
+++ b/field/devsim_field.py
@@ -10,16 +10,10 @@
import pickle
import ROOT
import numpy as np
-from scipy.interpolate import interp1d as p1d
-from scipy.interpolate import interp2d as p2d
-from scipy.interpolate import griddata
-from scipy.interpolate import LinearNDInterpolator as LNDI
-diff_res = 1e-5 # difference resolution in cm
+from util.math import *
-x_bin = 1000
-y_bin = 1000
-z_bin = 1000
+diff_res = 1e-5 # difference resolution in cm
class DevsimField:
def __init__(self, device_name, dimension, voltage, read_ele_num, l_z):
@@ -29,15 +23,39 @@ class DevsimField:
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)
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)
- self.set_potential(PotentialFile) #self.potential, self.x_efield, self.y_efield, self.z_efield
+ 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[]
+ def set_doping(self, DopingFile):
+ try:
+ with open(DopingFile,'rb') as file:
+ DopingNotUniform=pickle.load(file)
+ print("Doping file loaded for {}".format(self.name))
+ if DopingNotUniform['metadata']['dimension'] < self.dimension:
+ print("Doping dimension not match")
+ return
+ except FileNotFoundError:
+ print("Doping file not found, please run field simulation first")
+ print("or manually set the doping file")
+ return
+
+ if DopingNotUniform['metadata']['dimension'] == 1:
+ DopingUniform = get_common_interpolate_1d(DopingNotUniform)
+ elif DopingNotUniform['metadata']['dimension'] == 2:
+ DopingUniform = get_common_interpolate_2d(DopingNotUniform)
+ elif DopingNotUniform['metadata']['dimension'] == 3:
+ DopingUniform = get_common_interpolate_3d(DopingNotUniform)
+
+ self.Doping = DopingUniform
+
def set_potential(self, PotentialFile):
try:
with open(PotentialFile,'rb') as file:
@@ -54,9 +72,9 @@ class DevsimField:
if PotentialNotUniform['metadata']['dimension'] == 1:
PotentialUniform = get_common_interpolate_1d(PotentialNotUniform)
elif PotentialNotUniform['metadata']['dimension'] == 2:
- PotentialUniform =get_common_interpolate_2d(PotentialNotUniform)
+ PotentialUniform = get_common_interpolate_2d(PotentialNotUniform)
elif PotentialNotUniform['metadata']['dimension'] == 3:
- PotentialUniform =get_common_interpolate_3d(PotentialNotUniform)
+ PotentialUniform = get_common_interpolate_3d(PotentialNotUniform)
self.Potential = PotentialUniform
@@ -88,9 +106,9 @@ class DevsimField:
if TrappingRate_pNotUniform['metadata']['dimension'] == 1:
TrappingRate_pUniform = get_common_interpolate_1d(TrappingRate_pNotUniform)
elif TrappingRate_pNotUniform['metadata']['dimension'] == 2:
- TrappingRate_pUniform =get_common_interpolate_2d(TrappingRate_pNotUniform)
+ TrappingRate_pUniform = get_common_interpolate_2d(TrappingRate_pNotUniform)
elif TrappingRate_pNotUniform['metadata']['dimension'] == 3:
- TrappingRate_pUniform =get_common_interpolate_3d(TrappingRate_pNotUniform)
+ TrappingRate_pUniform = get_common_interpolate_3d(TrappingRate_pNotUniform)
self.TrappingRate_p = TrappingRate_pUniform
@@ -110,14 +128,27 @@ class DevsimField:
if TrappingRate_nNotUniform['metadata']['dimension'] == 1:
TrappingRate_nUniform = get_common_interpolate_1d(TrappingRate_nNotUniform)
elif TrappingRate_nNotUniform['metadata']['dimension'] == 2:
- TrappingRate_nUniform =get_common_interpolate_2d(TrappingRate_nNotUniform)
+ TrappingRate_nUniform = get_common_interpolate_2d(TrappingRate_nNotUniform)
elif TrappingRate_nNotUniform['metadata']['dimension'] == 3:
- TrappingRate_nUniform =get_common_interpolate_3d(TrappingRate_nNotUniform)
+ TrappingRate_nUniform = get_common_interpolate_3d(TrappingRate_nNotUniform)
self.TrappingRate_n = TrappingRate_nUniform
# DEVSIM dimension order: x, y, z
# RASER dimension order: z, x, y
+
+ def get_doping(self, x, y, z):
+ '''
+ input: position in um
+ output: doping in cm^-3
+ '''
+ x, y, z = x/1e4, y/1e4, z/1e4 # um to cm
+ if self.dimension == 1:
+ return self.Doping(z)
+ elif self.dimension == 2:
+ return self.Doping(z, x)
+ elif self.dimension == 3:
+ return self.Doping(z, x, y)
def get_potential(self, x, y, z):
'''
@@ -243,45 +274,6 @@ class DevsimField:
elif self.dimension == 3:
return self.TrappingRate_p(z, x, y)
-def get_common_interpolate_1d(data):
- values = data['values']
- points = data['points']
- f = p1d(points, values)
- return f
-
-def get_common_interpolate_2d(data):
- values = data['values']
- points_x = []
- points_y = []
- for point in data['points']:
- points_x.append(point[0])
- points_y.append(point[1])
- new_x = np.linspace(min(points_x), max(points_x), x_bin)
- new_y = np.linspace(min(points_y), max(points_y), y_bin)
- new_points = np.array(np.meshgrid(new_x, new_y)).T.reshape(-1, 2)
- new_values = griddata((points_x, points_y), values, new_points, method='linear')
- f = p2d(new_x, new_y, new_values)
- return f
-
-def get_common_interpolate_3d(data):
- values = data['values']
- points_x = []
- points_y = []
- points_z = []
- for point in data['points']:
- points_x.append(point[0])
- points_y.append(point[1])
- points_z.append(point[2])
- new_x = np.linspace(min(points_x), max(points_x), x_bin)
- new_y = np.linspace(min(points_y), max(points_y), y_bin)
- new_z = np.linspace(min(points_z), max(points_z), z_bin)
- new_points = np.array(np.meshgrid(new_x, new_y, new_z)).T.reshape(-1, 3)
- new_values = griddata((points_x, points_y, points_z), values, new_points, method='linear')
- lndi = LNDI(new_points, new_values)
- def f(x, y, z):
- point = [x, y, z]
- return lndi(point)
- return f
def linear_w_p(z, l_z):
if z >= l_z:
diff --git a/field/devsim_solve.py b/field/devsim_solve.py
index e2dd236b788b8ef6bab0caa9d0fee997624f8a68..bc648a6b23fbb576623b58d3aa05333a1ba52773 100644
--- a/field/devsim_solve.py
+++ b/field/devsim_solve.py
@@ -293,7 +293,11 @@ def milestone_save_1D(device, region, v, path):
metadata['voltage'] = v
metadata['dimension'] = 1
- for name in ['Potential', 'TrappingRate_p', 'TrappingRate_n']: # scalar field on mesh point (instead of on edge)
+ names = ['Potential', 'TrappingRate_p', 'TrappingRate_n']
+ 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:
data = {}
data['values'] = eval(name) # refer to the object with given name
@@ -327,7 +331,11 @@ def milestone_save_2D(device, region, v, path):
metadata['voltage'] = v
metadata['dimension'] = 2
- for name in ['Potential', 'TrappingRate_p', 'TrappingRate_n']: # scalar field on mesh point (instead of on edge)
+ names = ['Potential', 'TrappingRate_p', 'TrappingRate_n']
+ 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:
data = {}
data['values'] = eval(name) # refer to the object with given name
@@ -338,7 +346,6 @@ def milestone_save_2D(device, region, v, path):
pickle.dump(data, file)
-
def milestone_save_wf_2D(device, region, v, path,contact):
x = np.array(devsim.get_node_model_values(device=device, region=region, name="x")) # get x-node values
y = np.array(devsim.get_node_model_values(device=device, region=region, name="y")) # get y-node values
@@ -354,7 +361,6 @@ def milestone_save_wf_2D(device, region, v, path,contact):
draw2D(x,y,Potential,"Potential",v, path)
draw2D(x_mid,y_mid,ElectricField,"ElectricField",v, path)
-
dd = os.path.join(path, str(v),str(contact)+'V.dd')
devsim.write_devices(file=dd, type="tecplot")
@@ -382,7 +388,11 @@ def milestone_save_3D(device, region, v, path):
metadata['voltage'] = v
metadata['dimension'] = 3
- for name in ['Potential']: # scalar field on mesh point (instead of on edge)
+ names = ['Potential', 'TrappingRate_p', 'TrappingRate_n']
+ 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:
data = {}
data['values'] = eval(name) # refer to the object with given name
diff --git a/gen_signal/__init__.py b/gen_signal/__init__.py
index a4a63a362ced60c83869255528a7d4b387c7c172..af5d0aaf5ee383f6caf43689c730d8d7bd415560 100644
--- a/gen_signal/__init__.py
+++ b/gen_signal/__init__.py
@@ -1,204 +1,8 @@
-#!/usr/bin/env python3
-# -*- encoding: utf-8 -*-
-'''
-@Description: The main program of Raser induced current simulation
-@Date : 2024/02/20 18:12:26
-@Author : tanyuhang, Chenxi Fu
-@version : 2.0
-'''
-import sys
-import os
-import array
-import time
-import subprocess
-import ROOT
-
-from field import build_device as bdv
-from particle import g4simulation as g4s
-from field import devsim_field as devfield
-from current import cal_current as ccrt
-from elec import ele_readout as rdout
-from elec import ngspice_set_input as ngsip
-from elec import ngspice as ng
-
-from . import draw_save
-from util.output import output
-
-import json
-
-import random
-
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
- DevsimCal -- 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
- draw_plots -- Draw electric field, drift path and energy deposition
- Modify:
- ---------
- 2021/09/02
- """
- start = time.time()
-
- det_name = kwargs['det_name']
- my_d = bdv.Detector(det_name)
-
- if kwargs['voltage'] != None:
- voltage = float(kwargs['voltage'])
- else:
- voltage = float(my_d.voltage)
-
- if kwargs['absorber'] != None:
- absorber = kwargs['absorber']
- else:
- absorber = my_d.absorber
-
- if kwargs['amplifier'] != None:
- amplifier = kwargs['amplifier']
- else:
- amplifier = my_d.amplifier
-
- 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)
-
- if kwargs['scan'] != None:
- geant4_json = "./setting/absorber/" + absorber + ".json"
- with open(geant4_json) as f:
- g4_dic = json.load(f)
-
- total_events = int(g4_dic['total_events'])
- for i in range(kwargs['scan']):
- # TODO: change this into multithread
- instance_number = i
- g4_seed = instance_number * total_events
- my_g4p = g4s.Particles(my_d, absorber, g4_seed)
- batch_loop(my_d, my_f, my_g4p, amplifier, g4_seed, total_events, instance_number)
- del my_g4p
- return
-
- else:
- g4_seed = random.randint(0,1e7)
- my_g4p = g4s.Particles(my_d, absorber, g4_seed)
-
- 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)
+ label = kwargs['label']
- if 'ngspice' in amplifier:
- save_current(my_d, my_current, my_f = devfield.DevsimField(my_d.device, my_d.dimension, voltage, 1, my_d.l_z), key=None)
- input_p=ngsip.set_input(my_current, my_d, key=None)
- input_c=','.join(input_p)
- ng.ngspice_t0(input_c, input_p)
- subprocess.run(['ngspice -b -r t0.raw output/T0_tmp.cir'], shell=True)
- ng.plot_waveform()
+ if label == 'signal':
+ from . import gen_signal_main
+ gen_signal_main.main(kwargs)
else:
- ele_current = rdout.Amplifier(my_current, amplifier)
- draw_save.draw_plots(my_d,ele_current,my_f,my_g4p,my_current)
-
- del my_f
- end = time.time()
- print("total_time:%s"%(end-start))
-
-
-
-def batch_loop(my_d, my_f, my_g4p, amplifier, g4_seed, total_events, instance_number):
- """
- Description:
- Batch run some events to get time resolution
- Parameters:
- ---------
- start_n : int
- Start number of the event
- end_n : int
- end number of the event
- detection_efficiency: float
- The ration of hit particles/total_particles
- @Returns:
- ---------
- None
- @Modify:
- ---------
- 2021/09/07
- """
- path = output(__file__, my_d.det_name, 'batch')
- if "plugin" in my_d.det_model:
- draw_save.draw_ele_field(my_d,my_f,"xy",my_d.det_model,my_d.l_z*0.5,path)
- else:
- draw_save.draw_ele_field_1D(my_d,my_f,path)
- draw_save.draw_ele_field(my_d,my_f,"xz",my_d.det_model,my_d.l_y*0.5,path)
-
- start_n = instance_number * total_events
- end_n = (instance_number + 1) * total_events
-
- effective_number = 0
- for event in range(start_n,end_n):
- print("run events number:%s"%(event))
- if len(my_g4p.p_steps[event-start_n]) > 5:
- effective_number += 1
- my_current = ccrt.CalCurrentG4P(my_d, my_f, my_g4p, event-start_n)
- ele_current = rdout.Amplifier(my_current, amplifier)
- draw_save.save_signal_time_resolution(my_d,event,ele_current,my_g4p,start_n,my_f)
- del ele_current
- detection_efficiency = effective_number/(end_n-start_n)
- print("detection_efficiency=%s"%detection_efficiency)
-
-# TODO: change this to a method of CalCurrent
-def save_current(my_d,my_current,my_f,key):
- if key!=None:
- 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, )
- if os.path.exists(path):
- os.mkdir(path)
- L = eval("my_l.{}".format(key))
- #L is defined by different keys
- elif key==None:
- if "planar3D" in my_d.det_model or "planarRing" in my_d.det_model:
- path = os.path.join('output', 'PIN', my_d.det_name, )
- elif "lgad3D" in my_d.det_model:
- path = os.path.join('output', 'LGAD', my_d.det_name, )
- if not os.path.exists(path):
- os.makedirs(path)
-
- #L is defined by different keys
- time = array.array('d', [999.])
- current = array.array('d', [999.])
- fout = ROOT.TFile(os.path.join(path, "sim-current") + ".root", "RECREATE")
- t_out = ROOT.TTree("tree", "signal")
- t_out.Branch("time", time, "time/D")
- for i in range(my_current.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()
-
-
-
-
-
-
-
-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
+ raise NameError
diff --git a/gen_signal/gen_signal_main.py b/gen_signal/gen_signal_main.py
new file mode 100644
index 0000000000000000000000000000000000000000..1a7ea2821b7fc9dfd268b7a4a5096787535af62f
--- /dev/null
+++ b/gen_signal/gen_signal_main.py
@@ -0,0 +1,201 @@
+#!/usr/bin/env python3
+# -*- encoding: utf-8 -*-
+'''
+@Description: The main program of Raser induced current simulation
+@Date : 2024/02/20 18:12:26
+@Author : tanyuhang, Chenxi Fu
+@version : 2.0
+'''
+import sys
+import os
+import array
+import time
+import subprocess
+import ROOT
+
+from field import build_device as bdv
+from particle import g4simulation as g4s
+from field import devsim_field as devfield
+from current import cal_current as ccrt
+from elec import ele_readout as rdout
+from elec import ngspice_set_input as ngsip
+from elec import ngspice as ng
+
+from . import draw_save
+from util.output import output
+
+import json
+
+import random
+
+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
+ DevsimCal -- 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
+ draw_plots -- Draw electric field, drift path and energy deposition
+ Modify:
+ ---------
+ 2021/09/02
+ """
+ start = time.time()
+
+ det_name = kwargs['det_name']
+ my_d = bdv.Detector(det_name)
+
+ if kwargs['voltage'] != None:
+ voltage = float(kwargs['voltage'])
+ else:
+ voltage = float(my_d.voltage)
+
+ if kwargs['absorber'] != None:
+ absorber = kwargs['absorber']
+ else:
+ absorber = my_d.absorber
+
+ if kwargs['amplifier'] != None:
+ amplifier = kwargs['amplifier']
+ 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)
+
+ if kwargs['scan'] != None:
+ geant4_json = "./setting/absorber/" + absorber + ".json"
+ with open(geant4_json) as f:
+ g4_dic = json.load(f)
+
+ total_events = int(g4_dic['total_events'])
+ for i in range(kwargs['scan']):
+ # TODO: change this into multithread
+ instance_number = i
+ g4_seed = instance_number * total_events
+ my_g4p = g4s.Particles(my_d, absorber, g4_seed)
+ batch_loop(my_d, my_f, my_g4p, amplifier, g4_seed, total_events, instance_number)
+ del my_g4p
+ return
+
+ else:
+ g4_seed = random.randint(0,1e7)
+ my_g4p = g4s.Particles(my_d, absorber, g4_seed)
+
+ 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)
+
+ if 'ngspice' in amplifier:
+ save_current(my_d, my_current, my_f = devfield.DevsimField(my_d.device, my_d.dimension, voltage, 1, my_d.l_z), key=None)
+ input_p=ngsip.set_input(my_current, my_d, key=None)
+ input_c=','.join(input_p)
+ ng.ngspice_t0(input_c, input_p)
+ subprocess.run(['ngspice -b -r t0.raw output/T0_tmp.cir'], shell=True)
+ ng.plot_waveform()
+ else:
+ ele_current = rdout.Amplifier(my_current, amplifier)
+ draw_save.draw_plots(my_d,ele_current,my_f,my_g4p,my_current)
+
+ del my_f
+ end = time.time()
+ print("total_time:%s"%(end-start))
+
+
+
+def batch_loop(my_d, my_f, my_g4p, amplifier, g4_seed, total_events, instance_number):
+ """
+ Description:
+ Batch run some events to get time resolution
+ Parameters:
+ ---------
+ start_n : int
+ Start number of the event
+ end_n : int
+ end number of the event
+ detection_efficiency: float
+ The ration of hit particles/total_particles
+ @Returns:
+ ---------
+ None
+ @Modify:
+ ---------
+ 2021/09/07
+ """
+ path = output(__file__, my_d.det_name, 'batch')
+ if "plugin" in my_d.det_model:
+ draw_save.draw_ele_field(my_d,my_f,"xy",my_d.det_model,my_d.l_z*0.5,path)
+ else:
+ draw_save.draw_ele_field_1D(my_d,my_f,path)
+ draw_save.draw_ele_field(my_d,my_f,"xz",my_d.det_model,my_d.l_y*0.5,path)
+
+ start_n = instance_number * total_events
+ end_n = (instance_number + 1) * total_events
+
+ effective_number = 0
+ for event in range(start_n,end_n):
+ print("run events number:%s"%(event))
+ if len(my_g4p.p_steps[event-start_n]) > 5:
+ effective_number += 1
+ my_current = ccrt.CalCurrentG4P(my_d, my_f, my_g4p, event-start_n)
+ ele_current = rdout.Amplifier(my_current, amplifier)
+ draw_save.save_signal_time_resolution(my_d,event,ele_current,my_g4p,start_n,my_f)
+ del ele_current
+ detection_efficiency = effective_number/(end_n-start_n)
+ print("detection_efficiency=%s"%detection_efficiency)
+
+# TODO: change this to a method of CalCurrent
+def save_current(my_d,my_current,my_f,key):
+ if key!=None:
+ 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, )
+ if os.path.exists(path):
+ os.mkdir(path)
+ L = eval("my_l.{}".format(key))
+ #L is defined by different keys
+ elif key==None:
+ if "planar3D" in my_d.det_model or "planarRing" in my_d.det_model:
+ path = os.path.join('output', 'PIN', my_d.det_name, )
+ elif "lgad3D" in my_d.det_model:
+ path = os.path.join('output', 'LGAD', my_d.det_name, )
+ if not os.path.exists(path):
+ os.makedirs(path)
+
+ #L is defined by different keys
+ time = array.array('d', [999.])
+ current = array.array('d', [999.])
+ fout = ROOT.TFile(os.path.join(path, "sim-current") + ".root", "RECREATE")
+ t_out = ROOT.TTree("tree", "signal")
+ t_out.Branch("time", time, "time/D")
+ for i in range(my_current.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()
+
+
+
+
+
+
+
+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
diff --git a/motor/mdo_3034_control.py b/motor/mdo_3034_control.py
deleted file mode 100644
index 4edb4aaa36c32223e5e69b54ce20310ab922d30a..0000000000000000000000000000000000000000
--- a/motor/mdo_3034_control.py
+++ /dev/null
@@ -1,94 +0,0 @@
-import numpy as np
-import pandas as pd
-import pylab
-import re
-import time
-
-try:
- import visa
-except ImportError as err:
- print(err)
- exit()
-
-class MDO3034C:
- def __init__(self, resource_name):
- instlist=visa.ResourceManager()
- #print(instlist.list_resources())
- self.my_resource=instlist.open_resource(resource_name)
-
- def write_command(self,command):
- self.my_resource.write(command)
- time.sleep(0.5)
- def testIO(self):
- #self.my_resource.write('*RST')
- self.my_resource.write('*CLS')
- #self.my_resource.query_delay = 10.0
- message=self.my_resource.query('*IDN?')
- time.sleep(1)
- #print("ocsilloscope information:" + message)
- return message
-
- def readSet(self,ch,point_number):
- self.my_resource.write(":DATA:SOU "+ ch)
- self.my_resource.write(':DATA:START 1')
- self.my_resource.write(':DATA:STOP ' + point_number)
- self.my_resource.write(':WFMOutpre:ENCdg BINARY')
- self.my_resource.write(":WFMOutpre:BYT_Nr 1")
- self.my_resource.write(":HEADer 0")
- self.my_resource.write(":WFMOutpre?")
-
- def readOffset(self):
- self.write_command('WFMPRE:YMULT?')
- ymult = float(self.my_resource.read_raw())
- print("ymult = " + repr(ymult) + '\n')
- self.write_command('WFMPRE:YZERO?')
- yzero = float(self.my_resource.read_raw())
- print("yzero = " + repr(yzero) + '\n')
- self.write_command('WFMPRE:YOFF?')
- yoff = float(self.my_resource.read_raw())
- print("yoff = " + repr(yoff) + '\n')
- self.write_command('WFMPRE:XINCR?')
- xincr = float(self.my_resource.read_raw())
- print("xincr = " + repr(xincr) + '\n')
- self.write_command('WFMPRE:XZERO?')
- xzero = float(self.my_resource.read_raw())
- print("xzero = " + repr(xzero) + '\n')
- return ymult,yzero,yoff,xincr,xzero
-
- def readWave(self,ymult,yzero,yoff,xincr,xzero,point_num):
- #ymult,yzero,yoff,xincr,xzero=self.readOffset()
- self.my_resource.write('*CLS')
- self.my_resource.write("CURVE?")
- #file = open('test.txt','w+')
-
- ##########################################
- # receive data format:
- # #510000<data>\n receive 10000 data
- # #41000<data>\n receive 1000 data
- #########################################
- data = np.frombuffer(self.my_resource.read_raw(),dtype=np.int8,count=int(point_num),offset=len(str(point_num)) + 2)
- #np.savetxt('test.txt',data,fmt='%d',delimiter=',')
-
- print(data.size)
- Volts = (data - yoff) * ymult + yzero
- Time = np.arange(0, data.size, 1)
- Time = Time * xincr + xzero
-
- return Time, Volts
-
- def save_wave_data(self,time,voltage,filenames='./data.csv'):
- datadic = {'Time[ms]':time,'Voltage[mv]':voltage}
- dataform = pd.DataFrame(datadic,columns=['Time[ms]','Voltage[mv]'])
- dataform.to_csv(filenames,mode='a+',index=False)
-
-
- def plotWave(self, Time, Volts):
- pylab.plot(Time, Volts)
- pylab.show()
-
-def ReadInterface():
- rm = visa.ResourceManager()
- print(rm.list_resources())
- list_sources = rm.list_resources()
- return list_sources
-
diff --git a/motor/thread.py b/motor/thread.py
index 043a15fb93719cd04680fa125a4b5453588e8750..2b7a6b6a738dc8c137ac315fcd4c5a8257243292 100644
--- a/motor/thread.py
+++ b/motor/thread.py
@@ -1,5 +1,5 @@
import pymotor,stage,raser.motor.vitual_device as vitual_device
-import raser.motor.mdo_3034_control as MD
+import raser.control.mdo_3034_scope_control as MD
import time
from datetime import datetime
import numpy as np
diff --git a/particle/g4simulation.py b/particle/g4simulation.py
index c112ea7c53f802d9e0c7b4e43c26eab746fa3106..42396202909a4bd35b622b86d26c9f9e9a03e725 100644
--- a/particle/g4simulation.py
+++ b/particle/g4simulation.py
@@ -29,7 +29,7 @@ class Particles:
#use in pixel_detector
_randx = None
_randy = None
- def __init__(self, my_d, absorber, g4_seed):
+ def __init__(self, my_d, absorber, g4_seed = random.randint(0, 1e7)):
"""
Description:
Geant4 main process
diff --git a/root/__init__.py b/root/__init__.py
index 804037ae655831d228d6eaee957658bb5ea34668..fc45015038cf102ace020c9b54eb36863b56d110 100644
--- a/root/__init__.py
+++ b/root/__init__.py
@@ -1,60 +1,5 @@
-import sys
-import os
-import ROOT
-import csv
-
-def convert_csv_to_root(input_dir, output_dir, label):
- com_name = []
- for file in os.listdir(input_dir):
- if file.endswith('.csv'):
- com_name.append(file)
- for name in com_name:
- if label == 'sicar1.1.8' and not name.startswith('sicar1.1.8'):
- continue
- elif label == 'sicar1.1.8-1' and not name.startswith('sicar1.1.8-1_'):
- continue
- elif label == 'sicar1.1.8-2' and not name.startswith('sicar1.1.8-2_'):
- continue
-
- name = name.split('.csv')[0]
- input_file = os.path.join(input_dir, name + '.csv')
- output_file = os.path.join(output_dir, name + '.root')
-
- if name.endswith('iv'):
- if label=="itk_atlas18_data_v1":
- df = ROOT.RDF.MakeCsvDataFrame(input_file, True, '\t')
- elif label =="njupin_iv_v1":
- df = ROOT.RDF.MakeCsvDataFrame(input_file, True, ',')
- else:
- df = ROOT.RDF.MakeCsvDataFrame(input_file, True, ',')
- if label in ["itk_md8_data_v1","itk_atlas18_data_v1"]:
- df.Snapshot("myTree", output_file, {"Voltage_V", "Current_nA"})
- elif label in ['itk_atlas18_sim_v1','itk_md8_sim_v1']:
- df.Snapshot("myTree", output_file, {"Voltage", "Current"})
- elif label =="njupin_iv_v1":
- df.Snapshot("myTree", output_file, {"Current","Voltage"})
-
- else:
- df.Snapshot("myTree", output_file, {"Value","Reading"})
-
- if name.endswith('cv'):
- df = ROOT.RDF.MakeCsvDataFrame(input_file, True, ',')
- if label=="itk_md8_sim_v1":
- df.Snapshot("myTree", output_file, {"Voltage", "Capacitance"})
- elif label =="njupin_cv_v1":
- df.Snapshot("myTree", output_file, {"Voltage", "Capacitance"})
- else:
- df.Snapshot("myTree", output_file, {"Voltage", "Capacitance", "Capacitance^-2"})
-
- if name.endswith('Wfm'):
- df = ROOT.RDF.MakeCsvDataFrame(input_file, True, ',')
- df.Snapshot("myTree", output_file, {"Time", "Volt"})
-
-
- sys.stdout.write('Saved as {}\n'.format(output_file))
-
-
def main(kwargs):
+ from .convert import convert_csv_to_root
label = kwargs['label']
if label == 'sicar1.1.8':
diff --git a/root/convert.py b/root/convert.py
new file mode 100644
index 0000000000000000000000000000000000000000..09d92ada4b73d92f7c361b53cb70915fa7b8218c
--- /dev/null
+++ b/root/convert.py
@@ -0,0 +1,54 @@
+import sys
+import os
+import ROOT
+import csv
+
+def convert_csv_to_root(input_dir, output_dir, label):
+ com_name = []
+ for file in os.listdir(input_dir):
+ if file.endswith('.csv'):
+ com_name.append(file)
+ for name in com_name:
+ if label == 'sicar1.1.8' and not name.startswith('sicar1.1.8'):
+ continue
+ elif label == 'sicar1.1.8-1' and not name.startswith('sicar1.1.8-1_'):
+ continue
+ elif label == 'sicar1.1.8-2' and not name.startswith('sicar1.1.8-2_'):
+ continue
+
+ name = name.split('.csv')[0]
+ input_file = os.path.join(input_dir, name + '.csv')
+ output_file = os.path.join(output_dir, name + '.root')
+
+ if name.endswith('iv'):
+ if label=="itk_atlas18_data_v1":
+ df = ROOT.RDF.MakeCsvDataFrame(input_file, True, '\t')
+ elif label =="njupin_iv_v1":
+ df = ROOT.RDF.MakeCsvDataFrame(input_file, True, ',')
+ else:
+ df = ROOT.RDF.MakeCsvDataFrame(input_file, True, ',')
+ if label in ["itk_md8_data_v1","itk_atlas18_data_v1"]:
+ df.Snapshot("myTree", output_file, {"Voltage_V", "Current_nA"})
+ elif label in ['itk_atlas18_sim_v1','itk_md8_sim_v1']:
+ df.Snapshot("myTree", output_file, {"Voltage", "Current"})
+ elif label =="njupin_iv_v1":
+ df.Snapshot("myTree", output_file, {"Current","Voltage"})
+
+ else:
+ df.Snapshot("myTree", output_file, {"Value","Reading"})
+
+ if name.endswith('cv'):
+ df = ROOT.RDF.MakeCsvDataFrame(input_file, True, ',')
+ if label=="itk_md8_sim_v1":
+ df.Snapshot("myTree", output_file, {"Voltage", "Capacitance"})
+ elif label =="njupin_cv_v1":
+ df.Snapshot("myTree", output_file, {"Voltage", "Capacitance"})
+ else:
+ df.Snapshot("myTree", output_file, {"Voltage", "Capacitance", "Capacitance^-2"})
+
+ if name.endswith('Wfm'):
+ df = ROOT.RDF.MakeCsvDataFrame(input_file, True, ',')
+ df.Snapshot("myTree", output_file, {"Time", "Volt"})
+
+
+ sys.stdout.write('Saved as {}\n'.format(output_file))
diff --git a/spaceres/__init__.py b/spaceres/__init__.py
index f3c07cb0f9d34ac77831db38105086738de025a9..59b4b6e8a785215dbe8be42b51469446f7772dfa 100644
--- a/spaceres/__init__.py
+++ b/spaceres/__init__.py
@@ -1,17 +1,9 @@
import sys
import os
-#sys.path.insert(0, sys.path[0]+"/../")
-
-import ROOT
import time
import subprocess
-from . import telescope as tlcp
-#from . import test
-from particle.g4simulation import Particles
-from particle.geometry import R3dDetector
-from current.cal_current import CalCurrentPixel
-from util.output import create_path
+import ROOT
def main(kwargs):
label = kwargs['label']
@@ -19,103 +11,17 @@ def main(kwargs):
print("taichu_v1: ","first version of telescope simulation")
print("taichu_v2: ","second version of telescope simulation")
elif label.startswith("taichu_v1"):
- #paths = ['det_name=Taichu3', 'parfile=readjson/detector.json', 'geant4_model=pixel_detector', 'geant4_parfile=readjson/absorber.json', 'pixel_detector']
- dset = Setting() #read label.json instead of inputting path
- my_d = R3dDetector(dset) #remain the same
- my_f = 0
- my_g4p = Particles(my_d, dset) #remove my_f
- my_charge = CalCurrentPixel(my_d,my_f,my_g4p)
- if label.endswith("draw_charge"):
- draw_charge(my_charge)
- my_telescope = tlcp.telescope(my_d,my_charge)
- #tlcp.main(my_d)
+ from . import telescope_signal as tlcp
+ tlcp.main()
elif label.startswith("taichu_v2"):
- #virtual object
- class MyObject:
- pass
- #output
- res = []
- psize = []
- N = 25
- MaxSize = 25.
- for i in range(N):
- t_my_d = MyObject()
- t_my_d.seedcharge = 100
- t_my_d.p_x = MaxSize*(i+1)/N
- t_my_d.p_y = MaxSize*(i+1)/N
- t_my_d.p_z = 200.
- t_my_d.lt_z = [20000.,60000.,100000.,140000.,180000.,220000.]
- psize.append(t_my_d.p_x)
- res.append(tlcp.main(t_my_d))
-
- graph = ROOT.TGraph()
- for i in range(len(psize)):
- graph.SetPoint(i,psize[i],res[i])
-
- canvas = ROOT.TCanvas("canvas", "TGraph", 800, 600)
- graph.SetMarkerStyle(ROOT.kFullCircle)
- graph.GetYaxis().SetTitle("Resolution [um]")
- graph.GetYaxis().CenterTitle()
- graph.GetXaxis().SetTitle("Pixel Size [um]")
- graph.GetXaxis().CenterTitle()
-
- legend = ROOT.TLegend(0.27,0.67,0.62,0.80)
- legend.SetTextSize(0.04)
- legend.AddEntry(graph,label)
-
- canvas.SetGrid()
-
- graph.Draw("APL")
- legend.Draw()
-
- canvas.Draw()
- Name = "Res vs size"
- now = time.strftime("%Y_%m%d_%H%M")
- path = os.path.join("output/fig", str(now),'' )
- """ If the path does not exit, create the path"""
- if not os.access(path, os.F_OK):
- os.makedirs(path)
- canvas.SaveAs(path+Name+".png")
+ from . import telescope_signal as tlcp
+ tlcp.taichu_v2(label)
elif label.startswith("acts_v1"):
- python_script = "raser/spaceres/telescope_simulation.py"
-
- result = subprocess.run(["python3", python_script], stdout=subprocess.PIPE, stderr=subprocess.PIPE, text=True)
-
- if result.returncode == 0:
- print(f"Script '{python_script}' executed successfully.")
- else:
- print(f"Script '{python_script}' failed with return code:", result.returncode)
-
- stdout_output = result.stdout
- stderr_output = result.stderr
-
- print("Standard Output:")
- print(stdout_output)
-
- print("Standard Error:")
- print(stderr_output)
+ from . import telescope_acts as tlcp_acts
+ tlcp_acts.main()
+ elif label.startswith("g4"):
+ from . import telescope_g4
+ telescope_g4.main()
else:
raise NameError(label)
-def draw_charge(my_charge):
- path = os.path.join("output", "pixel",)
- create_path(path)
- c=ROOT.TCanvas("c","canvas1",1000,1000)
- c.cd()
- c.Update()
- c.SetLeftMargin(0.12)
- # c.SetTopMargin(0.12)
- c.SetRightMargin(0.12)
- c.SetBottomMargin(0.14)
- ROOT.gStyle.SetOptStat(ROOT.kFALSE)
- ROOT.gStyle.SetOptStat(0)
-
- my_charge.sum_charge.GetXaxis().SetNdivisions(510)
- my_charge.sum_charge.GetYaxis().SetNdivisions(505)
- my_charge.sum_charge.GetXaxis().SetTitle("X")
- my_charge.sum_charge.GetYaxis().SetTitle("Y")
-
- my_charge.sum_charge.Draw("lego")
- c.Update()
- c.SaveAs(path+"/Pixel_charge.pdf")
- c.SaveAs(path+"/Pixel_charge.root")
diff --git a/spaceres/telescope_simulation.py b/spaceres/telescope_acts.py
similarity index 93%
rename from spaceres/telescope_simulation.py
rename to spaceres/telescope_acts.py
index 80a2e6805dd4045e3a36d10da78fe575f81602bc..a007c578c63154f703cdfc079c9fd795ff4c9b32 100755
--- a/spaceres/telescope_simulation.py
+++ b/spaceres/telescope_acts.py
@@ -4,6 +4,9 @@ from pathlib import Path
import acts
import acts.examples
+
+from util.io_decorator import io_decorator as iod
+
from acts.examples.simulation import (
addParticleGun,
addDigitization,
@@ -30,9 +33,11 @@ from acts.examples.reconstruction import (
TrackSelectorConfig,
)
-u = acts.UnitConstants
+# get run result
+@iod
+def main():
+ u = acts.UnitConstants
-if "__main__" == __name__:
detector, trackingGeometry, decorators = acts.examples.TelescopeDetector.create(
positions=[20, 60, 100 ,140, 180, 220],
stereos = [0, 0, 0, 0, 0, 0],
@@ -42,10 +47,10 @@ if "__main__" == __name__:
field = acts.ConstantBField(acts.Vector3(0, 0, 1 * u.T))
- outputDir = Path.cwd() / "output/telescope_simulation"
+ outputDir = Path.cwd() / "output/spaceres/telescope_acts"
if not outputDir.exists():
outputDir.mkdir()
- inputDir = Path.cwd() / "setting"
+ inputDir = Path.cwd() / "setting/acts"
rnd = acts.examples.RandomNumbers(seed=42)
@@ -148,4 +153,7 @@ if "__main__" == __name__:
outputDirRoot=outputDir / postfix,
)
'''
- s.run()
\ No newline at end of file
+ s.run()
+
+if "__main__" == __name__:
+ main()
\ No newline at end of file
diff --git a/spaceres/telescope_jiaqi.py b/spaceres/telescope_g4.py
similarity index 94%
rename from spaceres/telescope_jiaqi.py
rename to spaceres/telescope_g4.py
index 1bdfd52bd9860fb599b9295bb16ca8f45c57135b..8d464f11b1dbfcdc29457cc24e2de55a0b769f84 100755
--- a/spaceres/telescope_jiaqi.py
+++ b/spaceres/telescope_g4.py
@@ -584,45 +584,50 @@ class B2ActionInitialization(G4VUserActionInitialization):
self.SetUserAction(B2EventAction())
self.SetUserAction(B2SteppingAction())
-# Detect interactive mode (if no arguments) and define UI session
-ui = None
-if len(sys.argv) == 1:
- ui = G4UIExecutive(len(sys.argv), sys.argv)
-
-# Optionally: choose a different Random engine...
-# G4Random.setTheEngine(MTwistEngine())
-
-# Construct the default run manager
-runManager = G4RunManagerFactory.CreateRunManager(G4RunManagerType.Serial)
-
-# Set mandatory initialization classes
-runManager.SetUserInitialization(B2aDetectorConstruction())
-
-physicsList = FTFP_BERT()
-physicsList.RegisterPhysics(G4StepLimiterPhysics())
-runManager.SetUserInitialization(physicsList)
-
-# Set user action classes
-runManager.SetUserInitialization(B2ActionInitialization())
-
-# Initialize visualization
-visManager = G4VisExecutive()
-# G4VisExecutive can take a verbosity argument - see /vis/verbose guidance.
-# visManager = G4VisExecutive("Quiet");
-visManager.Initialize()
-
-# Get the pointer to the User Interface manager
-UImanager = G4UImanager.GetUIpointer()
-
-# Process macro or start UI session
-if ui == None:
- # batch mode
- command = "/control/execute ./cfg/ "
- fileName = sys.argv[1]
- UImanager.ApplyCommand(command+fileName)
-else:
- # interactive mode
- UImanager.ApplyCommand("/control/execute ./param_file/g4macro/init_vistelescope.mac")
- if ui.IsGUI():
- UImanager.ApplyCommand("/control/execute ./param_file/g4macro/gui.mac")
- ui.SessionStart()
+
+def main(batchMac = None):
+ # Detect interactive mode (if no arguments) and define UI session
+ ui = None
+ if batchMac == None:
+ ui = G4UIExecutive(1, [__file__])
+
+ # Optionally: choose a different Random engine...
+ # G4Random.setTheEngine(MTwistEngine())
+
+ # Construct the default run manager
+ runManager = G4RunManagerFactory.CreateRunManager(G4RunManagerType.Serial)
+
+ # Set mandatory initialization classes
+ runManager.SetUserInitialization(B2aDetectorConstruction())
+
+ physicsList = FTFP_BERT()
+ physicsList.RegisterPhysics(G4StepLimiterPhysics())
+ runManager.SetUserInitialization(physicsList)
+
+ # Set user action classes
+ runManager.SetUserInitialization(B2ActionInitialization())
+
+ # Initialize visualization
+ visManager = G4VisExecutive()
+ # G4VisExecutive can take a verbosity argument - see /vis/verbose guidance.
+ # visManager = G4VisExecutive("Quiet");
+ visManager.Initialize()
+
+ # Get the pointer to the User Interface manager
+ UImanager = G4UImanager.GetUIpointer()
+
+ # Process macro or start UI session
+ if ui == None:
+ # batch mode
+ command = "/control/execute ./param_file/g4macro/"
+ fileName = batchMac
+ UImanager.ApplyCommand(command+fileName)
+ else:
+ # interactive mode
+ UImanager.ApplyCommand("/control/execute ./param_file/g4macro/init_vistelescope.mac")
+ if ui.IsGUI():
+ UImanager.ApplyCommand("/control/execute ./param_file/g4macro/gui.mac")
+ ui.SessionStart()
+
+if __name__ == "__main__":
+ main(sys.argv[1])
\ No newline at end of file
diff --git a/spaceres/telescope.py b/spaceres/telescope_signal.py
similarity index 84%
rename from spaceres/telescope.py
rename to spaceres/telescope_signal.py
index fbcc272202f7a5df50d83c252ad0c25219d08c89..2d2514377aa5187f8d098e7234dc08b53ac7c8ca 100644
--- a/spaceres/telescope.py
+++ b/spaceres/telescope_signal.py
@@ -6,13 +6,20 @@ Description:
@Author : yiminghu
@version : 1.0
'''
-import ROOT
+
import time
import os
+
+import ROOT
import acts
import numpy as np
-class telescope:
+from particle.g4simulation import Particles
+from field.build_device import Detector
+from current.cal_current_diffuse import CalCurrentPixel
+from util.output import create_path
+
+class Telescope:
def __init__(self,my_d,my_c):
"""
Description:
@@ -128,7 +135,7 @@ class telescope:
Name = "fit"+str(DUT)
now = time.strftime("%Y_%m%d_%H%M")
- path = os.path.join("fig", str(now),'' )
+ path = os.path.join("output/spaceres/taichu_v1", str(now),'' )
#print(path)
""" If the path does not exit, create the path"""
@@ -220,7 +227,7 @@ class telescope:
label.Draw()
now = time.strftime("%Y_%m%d_%H%M")
- path = os.path.join("fig", str(now),'' )
+ path = os.path.join("output/spaceres/taichu_v1", str(now),'' )
#print(path)
""" If the path does not exit, create the path"""
@@ -376,7 +383,7 @@ class island:
#interface to generate simple examples for debugging
-class Test:
+class PixelHitTest:
def __init__(self,my_d):
self.event = []
@@ -444,11 +451,89 @@ class Test:
#print(t_list)
return t_list
-def main(my_d):
- my_c = Test(my_d)
- tel = telescope(my_d,my_c)
- return tel.Resolution_Tol[2][0]
-
+def draw_charge(my_charge):
+ path = os.path.join("output", "pixel",)
+ create_path(path)
+ c=ROOT.TCanvas("c","canvas1",1000,1000)
+ c.cd()
+ c.Update()
+ c.SetLeftMargin(0.12)
+ # c.SetTopMargin(0.12)
+ c.SetRightMargin(0.12)
+ c.SetBottomMargin(0.14)
+ ROOT.gStyle.SetOptStat(ROOT.kFALSE)
+ ROOT.gStyle.SetOptStat(0)
+
+ my_charge.sum_charge.GetXaxis().SetNdivisions(510)
+ my_charge.sum_charge.GetYaxis().SetNdivisions(505)
+ my_charge.sum_charge.GetXaxis().SetTitle("X")
+ my_charge.sum_charge.GetYaxis().SetTitle("Y")
+
+ my_charge.sum_charge.Draw("lego")
+ c.Update()
+ c.SaveAs(path+"/Pixel_charge.pdf")
+ c.SaveAs(path+"/Pixel_charge.root")
+
+def main():
+ my_d = Detector("TAICHU3") #remain the same
+ my_f = 0
+ my_g4p = Particles(my_d, my_d.absorber) #remove my_f
+ my_hit_charge = CalCurrentPixel(my_d,my_f,my_g4p)
+ draw_charge(my_hit_charge)
+ my_telescope_charge = Telescope(my_d,my_hit_charge)
+ my_hit_test = PixelHitTest(my_d)
+ my_telescope_test = Telescope(my_d,my_hit_test)
+
+def taichu_v2(label=""):
+ #virtual object
+ class MyObject:
+ pass
+ #output
+ res = []
+ psize = []
+ N = 25
+ MaxSize = 25.
+ for i in range(N):
+ t_my_d = MyObject()
+ t_my_d.seedcharge = 100
+ t_my_d.p_x = MaxSize*(i+1)/N
+ t_my_d.p_y = MaxSize*(i+1)/N
+ t_my_d.p_z = 200.
+ t_my_d.lt_z = [20000.,60000.,100000.,140000.,180000.,220000.]
+ psize.append(t_my_d.p_x)
+ my_hit_test = PixelHitTest(t_my_d)
+ my_telescope_test = Telescope(t_my_d,my_hit_test)
+ res.append(my_telescope_test.Resolution_Tol[2][0])
+
+ graph = ROOT.TGraph()
+ for i in range(len(psize)):
+ graph.SetPoint(i,psize[i],res[i])
+
+ canvas = ROOT.TCanvas("canvas", "TGraph", 800, 600)
+ graph.SetMarkerStyle(ROOT.kFullCircle)
+ graph.GetYaxis().SetTitle("Resolution [um]")
+ graph.GetYaxis().CenterTitle()
+ graph.GetXaxis().SetTitle("Pixel Size [um]")
+ graph.GetXaxis().CenterTitle()
+
+ legend = ROOT.TLegend(0.27,0.67,0.62,0.80)
+ legend.SetTextSize(0.04)
+ legend.AddEntry(graph,label)
+
+ canvas.SetGrid()
+
+ graph.Draw("APL")
+ legend.Draw()
+
+ canvas.Draw()
+ Name = "Res vs size"
+ now = time.strftime("%Y_%m%d_%H%M")
+ path = os.path.join("output/spaceres/taichu_v2", str(now),'' )
+ """ If the path does not exit, create the path"""
+ if not os.access(path, os.F_OK):
+ os.makedirs(path)
+ canvas.SaveAs(path+Name+".png")
+
if __name__ == '__main__':
start = time.time()
main()
diff --git a/util/io_decorator.py b/util/io_decorator.py
new file mode 100644
index 0000000000000000000000000000000000000000..3dcc74dd612e5d4960a949d57c1f35529af0cce8
--- /dev/null
+++ b/util/io_decorator.py
@@ -0,0 +1,31 @@
+import sys
+import io
+from contextlib import redirect_stdout, redirect_stderr
+import importlib.util
+import inspect
+
+def io_decorator(func):
+ def wrapper(*args, **kwargs):
+ try:
+ # Capture the standard output and error
+ stdout_buffer = io.StringIO()
+ stderr_buffer = io.StringIO()
+
+ with redirect_stdout(stdout_buffer), redirect_stderr(stderr_buffer):
+ func(*args, **kwargs)
+ # Get the output
+ stdout_output = stdout_buffer.getvalue()
+ stderr_output = stderr_buffer.getvalue()
+
+ print(f"Function '{func.__name__}' executed successfully.")
+ print("Standard Output:")
+ print(stdout_output)
+
+ if stderr_output:
+ print("Standard Error:")
+ print(stderr_output)
+
+ except Exception as e:
+ print(f"Function '{func.__name__}' failed with an exception:", e)
+
+ return wrapper
diff --git a/util/math.py b/util/math.py
new file mode 100644
index 0000000000000000000000000000000000000000..acd31bf7b66f3df979bffd20ea7d7800736b58fb
--- /dev/null
+++ b/util/math.py
@@ -0,0 +1,97 @@
+#!/usr/bin/env python3
+# -*- encoding: utf-8 -*-
+
+'''
+Description:
+ Math Objects
+@Date : 2024/09/19 20:57:33
+@Author : Chenxi Fu
+@version : 1.0
+'''
+
+import math
+
+import numpy as np
+from scipy.interpolate import interp1d as p1d
+from scipy.interpolate import interp2d as p2d
+from scipy.interpolate import griddata
+from scipy.interpolate import LinearNDInterpolator as LNDI
+
+x_bin = 1000
+y_bin = 1000
+z_bin = 1000
+
+class Vector:
+ def __init__(self,a1,a2,a3):
+ self.components = [a1,a2,a3]
+
+ def cross(self,Vector_b):
+ """ Get vector cross product of self and another Vector"""
+ o1 = self.components[1]*Vector_b.components[2]-self.components[2]*Vector_b.components[1]
+ o2 = self.components[2]*Vector_b.components[0]-self.components[0]*Vector_b.components[2]
+ o3 = self.components[0]*Vector_b.components[1]-self.components[1]*Vector_b.components[0]
+ return Vector(o1,o2,o3)
+
+ def get_length(self):
+ " Return length of self"
+ return math.sqrt(self.components[0]*self.components[0]+self.components[1]*self.components[1]+self.components[2]*self.components[2])
+
+ def add(self,Vector_b):
+ " Return the sum of two Vectors. eg: [1,2,3]+[1,2,3] = [2,4,6]"
+ o1 = self.components[0]+Vector_b.components[0]
+ o2 = self.components[1]+Vector_b.components[1]
+ o3 = self.components[2]+Vector_b.components[2]
+ return Vector(o1,o2,o3)
+
+ def sub(self,Vector_b):
+ " Return the subtraction of two Vectors. eg: [1,2,3]-[1,2,3] = [0,0,0]"
+ o1 = self.components[0]-Vector_b.components[0]
+ o2 = self.components[1]-Vector_b.components[1]
+ o3 = self.components[2]-Vector_b.components[2]
+ return Vector(o1,o2,o3)
+
+ def mul(self,k):
+ " Return Vector multiplied by number. eg: 2 * [1,2,3] = [2,4,6]"
+ return Vector(self.components[0]*k,self.components[1]*k,self.components[2]*k)
+
+
+def get_common_interpolate_1d(data):
+ values = data['values']
+ points = data['points']
+ f = p1d(points, values)
+ return f
+
+def get_common_interpolate_2d(data):
+ values = data['values']
+ points_x = []
+ points_y = []
+ for point in data['points']:
+ points_x.append(point[0])
+ points_y.append(point[1])
+ new_x = np.linspace(min(points_x), max(points_x), x_bin)
+ new_y = np.linspace(min(points_y), max(points_y), y_bin)
+ new_points = np.array(np.meshgrid(new_x, new_y)).T.reshape(-1, 2)
+ new_values = griddata((points_x, points_y), values, new_points, method='linear')
+ f = p2d(new_x, new_y, new_values)
+ return f
+
+def get_common_interpolate_3d(data):
+ values = data['values']
+ points_x = []
+ points_y = []
+ points_z = []
+ for point in data['points']:
+ points_x.append(point[0])
+ points_y.append(point[1])
+ points_z.append(point[2])
+
+ new_x = np.linspace(min(points_x), max(points_x), x_bin)
+ new_y = np.linspace(min(points_y), max(points_y), y_bin)
+ new_z = np.linspace(min(points_z), max(points_z), z_bin)
+ new_points = np.array(np.meshgrid(new_x, new_y, new_z)).T.reshape(-1, 3)
+ new_values = griddata((points_x, points_y, points_z), values, new_points, method='linear')
+ lndi = LNDI(new_points, new_values)
+ def f(x, y, z):
+ point = [x, y, z]
+ return lndi(point)
+ return f
\ No newline at end of file