From 9b00600eeb7fe346a1c35b73741db867a6fb58c6 Mon Sep 17 00:00:00 2001
From: Yanpeng Li <lyp20001113@163.com>
Date: Fri, 27 Sep 2024 10:37:03 +0800
Subject: [PATCH] Change for CFLM, Version:2

---
 cflm/__init__.py          |  29 ++++++++
 cflm/cflm.json            |  62 ++++++++++++++++
 cflm/cflm_Volt_Curr.py    | 106 +++++++++++++++++++++++++++
 cflm/energy_deposition.py |  36 ++++++++++
 cflm/get_signal.py        |  87 ++++++++++++++++++++++
 cflm/ucsc.cir             | 147 ++++++++++++++++++++++++++++++++++++++
 6 files changed, 467 insertions(+)
 create mode 100644 cflm/__init__.py
 create mode 100644 cflm/cflm.json
 create mode 100644 cflm/cflm_Volt_Curr.py
 create mode 100644 cflm/energy_deposition.py
 create mode 100644 cflm/get_signal.py
 create mode 100644 cflm/ucsc.cir

diff --git a/cflm/__init__.py b/cflm/__init__.py
new file mode 100644
index 0000000..56030ac
--- /dev/null
+++ b/cflm/__init__.py
@@ -0,0 +1,29 @@
+import logging
+from . import cflm
+from . import RootPlot_1D
+from . import get_signal
+from . import cflm_Volt_Curr
+from . import energy_deposition
+
+def main(kwargs):
+    label = kwargs['label']
+    verbose = kwargs['verbose'] 
+
+    if verbose == 1: # -v 
+        logging.basicConfig(level=logging.INFO)
+    if verbose == 2: # -vv 
+        logging.basicConfig(level=logging.DEBUG)
+
+    logging.info('This is INFO messaage')
+    logging.debug('This is DEBUG messaage')
+
+    if label == 'cflm':
+       cflm.main()
+    if label == 'RootPlot_1D':
+       RootPlot_1D.RootPlot_1D()
+    if label == 'GetSignal':
+       get_signal.get_signal()
+    if label == 'GetVolCur':
+       cflm_Volt_Curr.getVolCur()
+    if label == 'GetEdep':
+       energy_deposition.getedep()
diff --git a/cflm/cflm.json b/cflm/cflm.json
new file mode 100644
index 0000000..e2224db
--- /dev/null
+++ b/cflm/cflm.json
@@ -0,0 +1,62 @@
+{
+    "object" : {
+        "elemental" : {
+            "pipe" : {
+                "name" : "pipe",
+                "material" : "G4_Cu",
+                "Rmin" : 28,
+                "Rmax" :31,
+                "Pipe_Z" : 200,
+                "PipeSphi" : 0,
+                "PipeDphi" : 180,
+                "colour" : [0, 1, 0],
+                "position_x" : 0,
+                "position_y" : 0,
+                "position_z" : 100
+    
+            }
+
+        },
+
+        "binary_compounds" : {
+            "detector" : {
+                "name" : "detector",
+                "material_1" : "Si",
+                "material_2" : "C",
+                "compound_name" :"SiC",
+                "density" : 3.2,
+                "natoms_1" : 50,
+                "natoms_2" : 50,
+                "side_x" : 50,                
+                "side_y" : 0.1,
+                "side_z" : 200,
+                "colour" : [0,0.5,0.8],
+                "position_x" : -31.05,
+                "position_y" : 0,
+                "position_z" : 100
+            }
+    }
+},
+
+     
+    "world"    : "G4_Galactic",
+    
+    "BeamOn"   : 1000,
+    
+    "NumofGun" : 1,
+    "par_type" : "e-",
+    "par_energy" : 24,
+    "par_in" : [-26.9, 0, -100],
+    "par_direct" : [[-0.0042, 0, 1]],
+    
+    "maxStep" : 0.5,
+
+    "DetModule" : "NJU-PIN.json",
+
+    "vis" : 0,
+
+    "EdepBaseName" : "110mm.root",
+    "PosBaseName"  : "test.txt",
+    "CurrentName"  : "Current_test.root" 
+
+}
diff --git a/cflm/cflm_Volt_Curr.py b/cflm/cflm_Volt_Curr.py
new file mode 100644
index 0000000..9f3b7b3
--- /dev/null
+++ b/cflm/cflm_Volt_Curr.py
@@ -0,0 +1,106 @@
+#!/usr/bin/env python3
+
+import ROOT
+import sys
+import numpy
+import json
+
+def read_file_voltage(file_path,file_name):
+    with open(file_path + '/' + file_name) as f:
+        lines = f.readlines()
+        time_v,volt = [],[]
+
+        for line in lines:
+            time_v.append(float(line.split()[0])*1e9)
+            volt.append(float(line.split()[1])*1e3)
+
+    time_v = numpy.array(time_v ,dtype='float64')
+    volt = numpy.array(volt,dtype='float64')
+
+    return time_v,volt
+
+def read_file_current(file_path,file_name):
+    with open(file_path + '/' + file_name) as f:
+        lines = f.readlines()
+        time_c,curr = [],[]
+
+        for line in lines:
+            time_c.append(float(line.split()[0])*1e9)
+            curr.append(float(line.split()[1])*1e6)
+
+    time_c = numpy.array(time_c ,dtype='float64')
+    curr = numpy.array(curr, dtype='float64')
+
+    return time_c, curr
+
+
+def getVolCur():
+    file_path = './raser/cflm/output'
+    
+    geant4_json = "./raser/cflm/cflm.json"
+    with open(geant4_json) as f:
+         g4_dic = json.load(f)
+
+    file_name_v = g4_dic['CurrentName'].split('.')[0] + '.raw'
+    file_name_c = g4_dic['CurrentName'].split('.')[0] + '_pwl_current.txt'
+    
+    com_name_v = file_name_v.split('.')[0]
+    com_name_c = file_name_c.split('.')[0]
+
+    time_v, volt, time_c, curr = [], [], [], []
+
+    time_v, volt = read_file_voltage(file_path,file_name_v)
+    length_v = len(time_v)
+    time_c, curr = read_file_current(file_path,file_name_c)
+    length_c = len(time_c)
+
+
+    ROOT.gROOT.SetBatch()
+        
+    c = ROOT.TCanvas('c','c',4000,2000)
+    
+    pad1 = ROOT.TPad("pad1", "pad1", 0.05, 0.05, 0.45, 0.95)
+    pad2 = ROOT.TPad("pad2", "pad2", 0.55, 0.05, 0.95, 0.95)
+
+    pad1.Draw()
+    pad2.Draw()
+    
+    pad1.cd()
+    f1 = ROOT.TGraph(length_c, time_c, curr)
+    f1.SetTitle(' ')
+    f1.SetLineColor(2)
+    f1.SetLineWidth(2)
+    f1.GetXaxis().SetTitle('Time [ns]')
+    f1.GetXaxis().SetLimits(0,10)
+    f1.GetXaxis().CenterTitle()
+    f1.GetXaxis().SetTitleSize(0.05)
+    f1.GetXaxis().SetTitleOffset(0.8)
+    f1.GetYaxis().SetTitle('Current [uA]')
+    f1.GetYaxis().SetLimits(0,-5)
+    f1.GetYaxis().CenterTitle()
+    f1.GetYaxis().SetTitleSize(0.07)
+    f1.GetYaxis().SetTitleOffset(0.7)
+    f1.Draw('AL')
+    pad1.Update()
+
+    pad2.cd()
+    f2 = ROOT.TGraph(length_v, time_v, volt)
+    f2.SetTitle(' ')
+    f2.SetLineColor(2)
+    f2.SetLineWidth(2)
+    f2.GetXaxis().SetTitle('Time [ns]')
+    f2.GetXaxis().SetLimits(0,10)
+    f2.GetXaxis().CenterTitle()
+    f2.GetXaxis().SetTitleSize(0.05)
+    f2.GetXaxis().SetTitleOffset(0.8)
+    f2.GetYaxis().SetTitle('Voltage [mV]')
+    f2.GetYaxis().SetLimits(0,-5)
+    f2.GetYaxis().CenterTitle()
+    f2.GetYaxis().SetTitleSize(0.07)
+    f2.GetYaxis().SetTitleOffset(0.7)
+    f2.Draw('AL')
+    pad2.Update()
+
+
+    c.SaveAs(f"Current_Voltage_{g4_dic['CurrentName'].split('.')[0]}.pdf")
+    
diff --git a/cflm/energy_deposition.py b/cflm/energy_deposition.py
new file mode 100644
index 0000000..20a03ae
--- /dev/null
+++ b/cflm/energy_deposition.py
@@ -0,0 +1,36 @@
+import ROOT
+import json
+
+def getedep():
+
+    up = 25
+    low = 0
+    nbins = 25
+
+    geant4_json = "./raser/cflm/cflm.json"
+    with open(geant4_json) as f:
+         g4_dic = json.load(f)
+
+    file = ROOT.TFile(f"raser/cflm/output/{g4_dic['EdepBaseName']}")
+    tree = file.Get("cflm")
+    branch = tree.GetBranch("Edetector")
+
+    data = []
+
+    for event in tree:
+        value = branch.GetLeaf("Edetector").GetValue()
+        data.append(value)
+
+    histogram = ROOT.TH1F("histogram", "Histogram of Energy deposition", nbins, low, up)
+
+    histogram.SetXTitle("Energy deposition(MeV)")
+    histogram.SetYTitle("Events/{:.3f}".format((up - low) / nbins))
+
+    for value in data:
+        histogram.Fill(value)
+
+    c = ROOT.TCanvas("c", "c", 800, 600)
+
+    histogram.Draw()
+    c.SaveAs(f"raser/cflm/output/{g4_dic['EdepBaseName'].split('.')[0]}.png")
+
diff --git a/cflm/get_signal.py b/cflm/get_signal.py
new file mode 100644
index 0000000..4cfd011
--- /dev/null
+++ b/cflm/get_signal.py
@@ -0,0 +1,87 @@
+#!/usr/bin/env python3
+# -*- encoding: utf-8 -*-
+
+import sys
+import os
+import array
+import time
+import subprocess
+import ROOT
+
+from field import build_device as bdv
+from . import cflm
+from field import devsim_field as devfield
+from current import cal_current as ccrt
+from elec.set_pwl_input import set_pwl_input as pwlin
+
+from util.output import output
+
+import json
+
+import re
+import numpy
+
+def get_signal():
+
+    geant4_json = "./raser/cflm/cflm.json"
+    with open(geant4_json) as f:
+         g4_dic = json.load(f)
+
+    detector_json = "./setting/detector/"
+    with open(os.path.join(detector_json , g4_dic['DetModule'])) as q:
+         det_dic = json.load(q)
+
+    start = time.time()
+
+    det_name = det_dic['det_name']
+    my_d = bdv.Detector(det_name)
+    voltage = det_dic['bias']['voltage']
+    amplifier = det_dic['amplifier']
+
+    print(my_d.device)
+    print(voltage)
+    
+    my_f = devfield.DevsimField(my_d.device, my_d.dimension, voltage, 1, my_d.l_z)
+
+    my_g4p = cflm.cflmG4Particles(my_d)
+
+    my_current = ccrt.CalCurrentG4P(my_d, my_f, my_g4p, 0)
+
+    if 'ngspice' in amplifier:
+        save_current(my_d, my_current, g4_dic, my_f = devfield.DevsimField(my_d.device, my_d.dimension, voltage, 1, my_d.l_z))
+        pwlin(f"raser/cflm/output/{g4_dic['CurrentName'].split('.')[0]}_pwl_current.txt", 'raser/cflm/ucsc.cir', 'raser/cflm/output/')
+        subprocess.run([f"ngspice -b -r ./xxx.raw raser/cflm/output/ucsc_tmp.cir"], shell=True)
+    
+    del my_f
+    end = time.time()
+    print("total_time:%s"%(end-start))
+    
+def save_current(my_d, my_current, g4_dic, my_f):
+
+    time = array.array('d', [999.])
+    current = array.array('d', [999.])
+    fout = ROOT.TFile(os.path.join("raser/cflm/output/", g4_dic['CurrentName'].split('.')[0])  + ".root", "RECREATE")
+    t_out = ROOT.TTree("tree", "signal")
+    t_out.Branch("time", time, "time/D")
+    for i in range(my_f.read_ele_num):
+        t_out.Branch("current"+str(i), current, "current"+str(i)+"/D")
+        for j in range(my_current.n_bin):
+            current[0]=my_current.sum_cu[i].GetBinContent(j)
+            time[0]=j*my_current.t_bin
+            t_out.Fill()
+    t_out.Write()
+    fout.Close()
+   
+    file = ROOT.TFile(os.path.join("raser/cflm/output/", g4_dic['CurrentName'].split('.')[0])  + ".root", "READ")
+    tree = file.Get("tree")
+
+    pwl_file = open(os.path.join("raser/cflm/output/", f"{g4_dic['CurrentName'].split('.')[0]}_pwl_current.txt"), "w")
+
+    for i in range(tree.GetEntries()):
+       tree.GetEntry(i)
+       time_pwl = tree.time
+       current_pwl = tree.current0
+       pwl_file.write(str(time_pwl) + " " + str(current_pwl) + "\n")
+    
+    pwl_file.close()
+    file.Close()
diff --git a/cflm/ucsc.cir b/cflm/ucsc.cir
new file mode 100644
index 0000000..4a3ab78
--- /dev/null
+++ b/cflm/ucsc.cir
@@ -0,0 +1,147 @@
+ucsc circuit 
+* name of the circuit
+
+* definition of BFR840L3RHESD, as a subcircuit
+.subckt BFR840L3RHESD 1 2 3
+*
+Rcx 15 1 1.57895
+Rbx 25 2 1.92983
+Rex 35 3 0.0800447
+*
+CBEPAR 22 33 1.9449E-013
+CBCPAR 22 11 3.44161E-014
+CCEPAR 11 33 2.24848E-013
+LB    22 20 3.04259E-010
+LC   11 10  2.88058E-010
+CBEPCK 20 30  1E-014
+CBCPCK 20 10  1.5502E-014
+CCEPCK 10 30  1E-014
+LBX    20 25 9.04631E-011
+LEX   30 35 3.71422E-011
+LCX   10 15  9.15043E-011
+*
+R_CS_npn 55 5 500
+*
+D1 33 25 M_D1
+D2 5 25  M_D2
+*
+R_NBL_fdb 22 25 3.2
+R_PS 33 5 0.03
+RSUB 30 5 0.03
+*
+D3 5 15 M_D3
+D4 23 33 M_D4
+D5 23 15 M_D5
+*
+R_NBL_e11g 15 11 1.8
+*
+Q1 11 22 33 55 M_BFR840L3RHESD
+*
+.MODEL M_D1 D(
++ IS=3E-015
++ N=1
++ RS=2.846
++ CJO=4E-014)
+*
+.MODEL M_D2 D(
++ IS=3E-015
++ N=1
++ RS=4170
++ CJO=4.5E-014)
+*
+.MODEL M_D3 D(
++ IS=6.911E-016
++ N=1.1
++ RS=1380
++ CJO =9.5E-014)
+*
+.MODEL M_D4 D(
++ IS=3.5E-015
++ N=1
++ RS=0.2
++ CJO =3E-014)
+*
+.MODEL M_D5 D(
++ IS=3.5E-015
++ N=1.02
++ RS=4.7
++ CJO =3E-014)
+*
+.MODEL 	M_BFR840L3RHESD	NPN(
++	TNOM = 25
++	IS	=	2.429E-016
++	BF	=	765.7
++	NF	=	1.012
++	VAF	=	375.1
++	IKF	=	0.0819
++	ISE	=	8.827E-014
++	NE	=	2.8
++	BR	=	194
++	NR	=	0.998
++	VAR	=	1.596
++	IKR	=	0.015
++	ISC	=	1.165E-015
++	NC	=	2
++	RB	=	7.53378
++	IRB	=	0
++	RBM	=	2.1
++	RE	=	0.4405
++	RC	=	7.246
++	XTB	=	-2.276
++	EG	=	1.11
++	XTI	=	-1.233
++	CJE	=	2.23E-014
++	VJE	=	0.9214
++	MJE	=	0.5
++	TF	=	1.1E-012
++	XTF	=	5.582
++	VTF	=	0.6828
++	ITF	=	0.4491
++	PTF	=	0.0214
++	CJC	=	6.6E-015
++	VJC	=	0.7723
++	MJC	=	1.005
++	XCJC	=	0.4894
++	TR	=	1E-010
++	CJS	=	1.147E-013
++	MJS	=	1.108
++   VJS =   0.6112
++	FC	=	0.578
++	KF	=	1.65E-011
++	AF	=	1.53)
+.ends BFR840L3RHESD
+
+I1 2 0 pulse(0 -10u 0 0.1n 1n 0.00000001n 20n 0)
+* input current source
+VCC 6 0 dc 2.25
+* VCC, DC source
+Rin 2 0 1MEG
+* resistance of the sensor
+C6 2 0 20p
+* capacitance of the sensor
+C1 2 3 3.3n
+x1 5 3 0 BFR840L3RHESD
+* the amplifier BFR840L3RHESD
+R1 3 4 475
+* feedback resistance
+R2 4 5 3k
+C2 4 5 3.3n
+C5 5 out 3.3n
+R4 out 0 50
+L1 5 6 47U
+R3 6 7 63.7
+C7 8 0 10n
+C8 7 0 10n
+C3 7 0 1n
+C4 7 0 1n
+L2 8 7 2.2u
+* some other devices
+
+.control
+tran 0.1p 20n
+* <step> <stopping time>
+wrdata ./raser/cflm/output/Current_test.raw v(out)
+* save raw data
+.endc
+
+.end
-- 
GitLab