diff --git a/bmos/__init__.py b/bmos/__init__.py
index 86844fca7a173a7abd077b8797d087621961e814..d83c957c18a435cd9d1488a3ecf6c79604c85a94 100644
--- a/bmos/__init__.py
+++ b/bmos/__init__.py
@@ -16,4 +16,15 @@ def main(kwargs):
if label == 'GetSignal':
from . import get_signal
get_signal.get_signal()
+
+ if label == 'histogram_signal':
+ from . import histogram_signal
+ histogram_signal.get_signal()
+ if label == 'one_histogram':
+ from . import histogram
+ histogram.main("one")
+
+ if label == 'all_histogram':
+ from . import histogram
+ histogram.main("all")
diff --git a/bmos/histogram.py b/bmos/histogram.py
new file mode 100644
index 0000000000000000000000000000000000000000..30b6c14066a43d233389ac8c5f0afe48ebd8c108
--- /dev/null
+++ b/bmos/histogram.py
@@ -0,0 +1,144 @@
+#!/usr/bin/env python3
+# -*- encoding: utf-8 -*-
+
+import sys
+import os
+import array
+import time
+import subprocess
+import re
+import json
+
+import ROOT
+ROOT.gROOT.SetBatch(True)
+import numpy
+
+from gen_signal import build_device as bdv
+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
+from . import bmos
+import numpy as np
+
+def mkdir(folder_name):
+ try:
+ os.makedirs(folder_name)
+ except Exception as e:
+ pass
+
+def readfile(path, root_name):
+ file = ROOT.TFile(os.path.join(path, root_name), "READ")
+ tree = file.Get("tree")
+
+ amplitudes = []
+ for i in range(tree.GetEntries()):
+ tree.GetEntry(i)
+ amplitudes.append(tree.amplitude)
+
+ file.Close()
+ return amplitudes
+
+def data(root_path, root_name):
+ file_name = os.path.splitext(root_name)[0]
+ amplitudes = readfile(root_path, root_name)
+ mean = np.mean(amplitudes)
+ std = np.std(amplitudes)
+ # amplitudes = [amplitude*17.925/0.138*0.8847486018957982 for amplitude in amplitudes if amplitude > mean - 3*std and amplitude < mean + 3*std]
+ amplitudes = [amplitude*17.925/0.138 for amplitude in amplitudes if amplitude > mean - 3*std and amplitude < mean + 3*std]
+
+ return amplitudes, file_name
+
+
+def draw(amplitudes, file_name):
+ c = ROOT.TCanvas( 'c', 'c', 8000, 6000 )
+ c.cd()
+
+ minsignal = float(min(amplitudes))
+ maxsignal = float(max(amplitudes))
+ binnum = 50
+ binwidth = (maxsignal - minsignal)/binnum
+ wave_graph = ROOT.TH1F('','', binnum + 2, minsignal - binwidth, maxsignal + binwidth)
+ for amplitude in amplitudes:
+ wave_graph.Fill(amplitude)
+
+ f1 = ROOT.TF1("f1", "landau")
+ wave_graph.Fit(f1)
+
+ wave_graph.SetTitle(file_name)
+ wave_graph.SetLineColor(1)
+ wave_graph.SetLineWidth(2)
+
+ wave_graph.GetXaxis().SetTitle('mV')
+ wave_graph.GetXaxis().SetTitleSize(0.05)
+ wave_graph.GetXaxis().SetTitleOffset(0.9)
+ wave_graph.GetYaxis().SetTitle('events')
+ wave_graph.GetYaxis().SetTitleSize(0.05)
+ wave_graph.GetYaxis().SetTitleOffset(0.9)
+ # wave_graph.GetXaxis().SetRangeUser(0, 1300)
+ # wave_graph.GetXaxis().SetLimits(0, 1300)
+ wave_graph.Draw()
+ maxevent = wave_graph.GetMaximum()
+
+ latex = ROOT.TLatex(0.4*minsignal + 0.6*maxsignal, 0.7*maxevent, f"MPV:{round(f1.GetParameter(1), 3)}mV")
+ latex.SetTextSize(0.05)
+ latex.SetTextColor(1)
+ latex.SetTextFont(42)
+ latex.Draw()
+
+ mkdir(os.path.join("/scratchfs/atlas/xiekaibo/sicar/raser/bmos/output/", 'histogram', 'pdf'))
+
+ c.SaveAs(os.path.join("/scratchfs/atlas/xiekaibo/sicar/raser/bmos/output/", 'histogram', 'pdf', f"{file_name}.pdf"))
+
+
+
+def main(choose):
+ tag = 'proton_80MeV'
+ # tag = ''
+ root_path = os.path.join("/scratchfs/atlas/xiekaibo/sicar/raser/bmos/output/histogram/root", tag)
+ root_names = os.listdir(root_path)
+
+ if choose == "all":
+ for root_name in root_names:
+ amplitudes, file_name = data(root_path, root_name)
+ draw(amplitudes, file_name)
+ elif choose == "one":
+ c = ROOT.TCanvas( 'c', 'c', 8000, 6000 )
+ his = []
+
+ for i in range(len(root_names)):
+ c.cd()
+
+ amplitudes, file_name = data(root_path, root_names[i])
+ # minsignal = float(min(amplitudes))
+ # maxsignal = float(max(amplitudes))
+ # binnum = 50
+ # binwidth = (maxsignal - minsignal)/binnum
+ # his.append(ROOT.TH1F('','', binnum + 2, minsignal - binwidth, maxsignal + binwidth))
+ his.append(ROOT.TH1F('','', 150, 0, 1500))
+
+ for amplitude in amplitudes:
+ his[i].Fill(amplitude)
+
+ his[i].Draw("same")
+ his[i].Fit("landau")
+
+ his[i].SetTitle(tag)
+ his[i].SetLineColor(1)
+ his[i].SetLineWidth(2)
+
+ his[i].GetXaxis().SetTitle('mV')
+ his[i].GetXaxis().SetTitleSize(0.05)
+ his[i].GetXaxis().SetTitleOffset(0.9)
+ # his[i].GetXaxis().SetRangeUser(0, 1300)
+ # his[i].GetXaxis().SetLimits(0, 1300)
+ his[i].GetYaxis().SetTitle('events')
+ his[i].GetYaxis().SetTitleSize(0.05)
+ his[i].GetYaxis().SetTitleOffset(0.9)
+ his[i].GetYaxis().SetRangeUser(0, 150)
+
+ mkdir(os.path.join("/scratchfs/atlas/xiekaibo/sicar/raser/bmos/output/", 'histogram', 'pdf', tag))
+ c.SaveAs(os.path.join("/scratchfs/atlas/xiekaibo/sicar/raser/bmos/output/", 'histogram', 'pdf', tag, f"{file_name}_all.pdf"))
+
+
+
\ No newline at end of file
diff --git a/bmos/histogram_signal.py b/bmos/histogram_signal.py
new file mode 100644
index 0000000000000000000000000000000000000000..5a222f591145df68e3003e3ab23c2b87438780b3
--- /dev/null
+++ b/bmos/histogram_signal.py
@@ -0,0 +1,171 @@
+#!/usr/bin/env python3
+# -*- encoding: utf-8 -*-
+
+import sys
+import os
+import array
+import time
+import subprocess
+import re
+import json
+
+import ROOT
+ROOT.gROOT.SetBatch(True)
+import numpy
+
+from gen_signal import build_device as bdv
+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
+from . import bmos
+
+def mkdir(folder_name):
+ try:
+ os.makedirs(folder_name)
+ except Exception as e:
+ pass
+
+def get_signal():
+ signal = []
+
+ geant4_json = "./setting/absorber/bmos.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']
+
+ my_f = devfield.DevsimField(my_d.device, my_d.dimension, voltage, my_d.read_out_contact, my_d.irradiation_flux)
+
+ my_g4p = bmos.bmosG4Particles(my_d)
+
+ output_path = "raser/bmos/output/"
+ tag = f"{g4_dic['par_type']}_{g4_dic['par_energy']}MeV_{g4_dic['par_num']}particle"
+ dirname = f"{g4_dic['par_type']}_{g4_dic['par_energy']}MeV"
+ root_name = f"{g4_dic['CurrentName'].split('.')[0]}_{tag}.root"
+ pwl_name = f"pwl{g4_dic['CurrentName'].split('.')[0]}_{tag}.txt"
+ filename_after_ngspice = f"UCSC_output_{tag}.raw"
+
+ for i in range(len(my_g4p.p_steps_current)):
+ my_current = ccrt.CalCurrentG4P(my_d, my_f, my_g4p, i)
+
+ save_current(my_current, output_path, root_name, pwl_name, 1)
+
+ pwlin(os.path.join(output_path, pwl_name), amplifier, output_path, filename_after_ngspice)
+ subprocess.run([f"ngspice -b -r ./xxx.raw raser/bmos/output/ucsc_tmp.cir"], shell=True)
+ time_v, volt = read_file_voltage(output_path, filename_after_ngspice)
+ signal.append(max(volt))
+
+ print(signal)
+ draw(output_path, signal, tag, dirname)
+
+ # draw(output_path, pwl_name, filename_after_ngspice, tag, totalengry)
+
+ end = time.time()
+ print("total_time:%s"%(end - start))
+
+def save_current(my_current, output_path, root_name, pwl_name, read_ele_num):
+ time = array.array('d', [999.])
+ current = array.array('d', [999.])
+
+ fout = ROOT.TFile(os.path.join(output_path, root_name), "RECREATE")
+ t_out = ROOT.TTree("tree", "signal")
+ t_out.Branch("time", time, "time/D")
+ for i in range(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(output_path, root_name), "READ")
+ tree = file.Get("tree")
+
+ pwl_file = open(os.path.join(output_path, pwl_name), "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()
+
+def read_file_voltage(file_path, file_name):
+ with open(os.path.join(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 draw(output_path, signal, tag, dirname):
+ ROOT.gROOT.SetBatch(True)
+ c = ROOT.TCanvas( 'c', 'c', 8000, 6000 )
+ # c.SetCanvasSize(800, 600)
+ # c.SetWindowSize(800, 600)
+ c.cd()
+ minsignal = float(min(signal))
+ maxsignal = float(max(signal))
+ binnum = 50
+ binwidth = (maxsignal - minsignal)/binnum
+ wave_graph = ROOT.TH1F('','', binnum + 2, minsignal - binwidth, maxsignal + binwidth)
+
+ mkdir(os.path.join(output_path, 'histogram', 'root', dirname))
+ file = ROOT.TFile(os.path.join(output_path, 'histogram', 'root', dirname, f"Histogram_{tag}.root"), "RECREATE")
+ amplitude = ROOT.TTree("tree", "amplitude")
+ amp =array.array('d', [0.0])
+ amplitude.Branch("amplitude", dirname, "amplitude/D")
+ # amplitude.Write()
+ # file.Close()
+
+ for sig in signal:
+ amp[0] = sig
+ amplitude.Fill()
+ wave_graph.Fill(sig)
+
+ # if ifnoise(time, volt, noisestd):
+ # num_noise += 1
+ # continue
+ # num_signal += 1
+ wave_graph.Draw()
+ wave_graph.SetTitle(' ')
+ wave_graph.SetLineColor(1)
+ wave_graph.SetLineWidth(2)
+
+ wave_graph.GetXaxis().SetTitle('mV')
+ wave_graph.GetXaxis().SetTitleSize(0.05)
+ wave_graph.GetXaxis().SetTitleOffset(0.9)
+ wave_graph.GetYaxis().SetTitle('events')
+ wave_graph.GetYaxis().SetTitleSize(0.05)
+ wave_graph.GetYaxis().SetTitleOffset(0.9)
+ # wave_graph.Fit("gaus")
+ wave_graph.Fit("landau")
+
+ # c.SaveAs(os.path.join(output_path, f"Histogram_{tag}.root"))
+ c.SaveAs(os.path.join(output_path, 'histogram', 'pdf', f"Histogram_{tag}.pdf"))
+
+ amplitude.Write()
+ wave_graph.Write()
+ file.Close()
+
+