diff --git a/__main__.py b/__main__.py
index c428ce797b848c0d03d0e2551d45b67679dc45b6..38a75e8cf62b219634791ca82a61b8edfed8c5ad 100755
--- a/__main__.py
+++ b/__main__.py
@@ -72,6 +72,10 @@ parser_tct.add_argument('-vol', '--voltage', type=str, help='bias voltage')
 parser_tct.add_argument('-amp', '--amplifier', type=str, help='amplifier')
 parser_tct.add_argument('-s', '--scan', type=int, help='instance number for scan mode')
 
+parser_bmos = subparsers.add_parser('bmos', help='Beam Monitor Online System')
+parser_bmos.add_argument('label', help='LABEL to identify BMOS simulations')
+parser_bmos.add_argument('-v', '--verbose', help='VERBOSE level', 
+                          action='count', default=0)
 
 args = parser.parse_args()
 
@@ -81,7 +85,7 @@ if len(sys.argv) == 1:
 
 kwargs = vars(args)
 
-submodules = ['asic', 'cce', 'cflm', 'current', 'draw', 'elec', 'field', 'fpga', 'gen_signal', 'particle', 'spaceres', 'tct', 'timeres']
+submodules = ['asic', 'cce', 'cflm', 'current', 'draw', 'elec', 'field', 'fpga', 'gen_signal', 'particle', 'spaceres', 'tct', 'timeres', 'bmos']
 
 submodule = kwargs['subparser_name']
 if submodule not in submodules:
diff --git a/bmos/__init__.py b/bmos/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..86844fca7a173a7abd077b8797d087621961e814
--- /dev/null
+++ b/bmos/__init__.py
@@ -0,0 +1,19 @@
+import logging
+
+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 == 'GetSignal':
+        from . import get_signal
+        get_signal.get_signal()
+
diff --git a/bmos/bmos.py b/bmos/bmos.py
new file mode 100644
index 0000000000000000000000000000000000000000..d81a67fd0bcffc68298430e2258dc0257d89d227
--- /dev/null
+++ b/bmos/bmos.py
@@ -0,0 +1,331 @@
+#!/usr/bin/env python3 
+import os
+import math
+
+import geant4_pybind as g4b
+import json
+import numpy as np
+import time
+
+#G4AnalysisManager = g4b.G4RootAnalysisManager
+
+X_position = []
+Z_position = []
+Y_position = []
+Particle = []
+
+
+class bmosG4Particles:
+
+    def __init__(self, my_d):
+
+        global s_eventIDs, s_edep_devices, s_p_steps, s_energy_steps
+        s_eventIDs, s_edep_devices, s_p_steps, s_energy_steps = [], [], [], []
+
+        geant4_json = "./setting/absorber/bmos.json"
+        with open(geant4_json) as f:
+             g4_dic = json.load(f)
+
+        self.geant4_model = g4_dic["geant4_model"]
+
+        runManager = g4b.G4RunManagerFactory.CreateRunManager(g4b.G4RunManagerType.Serial)
+        UImanager = g4b.G4UImanager.GetUIpointer()
+
+        physicsList = g4b.FTFP_BERT()
+        physicsList.SetVerboseLevel(1)
+        physicsList.RegisterPhysics(g4b.G4StepLimiterPhysics())
+        runManager.SetUserInitialization(physicsList)
+
+        detConstruction = DetectorConstruction(g4_dic)
+        runManager.SetUserInitialization(detConstruction)
+
+        actionInitialization = ActionInitialization(detConstruction,
+                                                    g4_dic['par_in'],
+                                                    g4_dic['par_direction'],
+                                                    g4_dic['par_type'],
+                                                    g4_dic['par_energy'],
+                                                    g4_dic['par_num'],
+                                                    )
+
+        runManager.SetUserInitialization(actionInitialization)
+
+        UImanager = g4b.G4UImanager.GetUIpointer()
+        UImanager.ApplyCommand('/run/initialize')
+
+        runManager.BeamOn(int(g4_dic['BeamOn']))
+
+        self.p_steps = s_p_steps
+        self.init_tz_device = 0
+        self.p_steps_current = [[[single_step[0]+my_d.l_x/2,
+                                  single_step[1]+my_d.l_y/2,
+                                  single_step[2]-self.init_tz_device]\
+                for single_step in p_step] for p_step in self.p_steps]
+
+        self.energy_steps = s_energy_steps
+        self.edep_devices = s_edep_devices
+
+        print("sum edep = ", sum(s_energy_steps[0]))
+        # time.sleep(10)
+
+        del s_eventIDs,s_edep_devices,s_p_steps,s_energy_steps
+ 
+
+class DetectorConstruction(g4b.G4VUserDetectorConstruction):
+
+    def __init__(self,g4_dic):
+        g4b.G4VUserDetectorConstruction.__init__(self)
+        self.solid = {}
+        self.logical = {}
+        self.physical = {}
+        self.checkOverlaps = True
+        self.create_world(g4_dic['world_type'], g4_dic['world_size'])
+
+        self.maxStep = g4_dic['maxStep']*g4b.um
+
+        # self.rotation = g4b.G4RotationMatrix()
+        # self.rotation.rotateZ(3*math.pi/2)
+
+        for object_type in g4_dic['object']:
+            if(object_type=="elemental"):
+                for every_object in g4_dic['object'][object_type]:
+                    self.create_elemental(g4_dic['object'][object_type][every_object])
+            if(object_type=="binary_compounds"):
+                for every_object in g4_dic['object'][object_type]:
+                    self.create_binary_compounds(g4_dic['object'][object_type][every_object])
+
+        self.fStepLimit = g4b.G4UserLimits(self.maxStep)
+        self.logical["detector"].SetUserLimits(self.fStepLimit)
+
+    def create_world(self, world_type, world_size):
+
+        self.nist = g4b.G4NistManager.Instance()
+        material = self.nist.FindOrBuildMaterial(world_type)  
+        self.solid['world'] = g4b.G4Box("world",
+                                        world_size*g4b.um,
+                                        world_size*g4b.um,
+                                        world_size*g4b.um)
+        self.logical['world'] = g4b.G4LogicalVolume(self.solid['world'], 
+                                                    material, 
+                                                    "world")
+        self.physical['world'] = g4b.G4PVPlacement(None, 
+                                                   g4b.G4ThreeVector(0,0,0), 
+                                                   self.logical['world'], 
+                                                   "world", None, False, 
+                                                   0,self.checkOverlaps)
+        
+        self.logical['world'].SetVisAttributes(g4b.G4VisAttributes.GetInvisible())
+
+    def create_elemental(self,object): 
+        name = object['name']
+        material_type = self.nist.FindOrBuildMaterial(object['material'],
+                                                    False)
+        translation = g4b.G4ThreeVector(object['position_x']*g4b.um, object['position_y']*g4b.um, object['position_z']*g4b.um)
+        visual = g4b.G4VisAttributes(g4b.G4Color(object['colour'][0],object['colour'][1],object['colour'][2]))
+        mother = self.physical['world']
+        sidex = object['side_x']*g4b.um
+        sidey = object['side_y']*g4b.um
+        sidez = object['side_z']*g4b.um
+        self.solid[name] = g4b.G4Box(name, sidex/2., sidey/2., sidez/2.)
+        
+        self.logical[name] = g4b.G4LogicalVolume(self.solid[name], 
+                                                material_type, 
+                                                name)
+        self.physical[name] = g4b.G4PVPlacement(None,translation,                                                
+                                                name,self.logical[name],
+                                                mother, False, 
+                                                0,self.checkOverlaps)
+        self.logical[name].SetVisAttributes(visual)     
+
+    def create_binary_compounds(self,object):
+        name = object['name']
+        material_1 = self.nist.FindOrBuildElement(object['material_1'],False)
+        material_2 = self.nist.FindOrBuildElement(object['material_2'],False)
+        material_density = object['density']*g4b.g/g4b.cm3
+        compound=g4b.G4Material(object['compound_name'],material_density,2) 
+        compound.AddElement(material_1,object['natoms_1']*g4b.perCent)
+        compound.AddElement(material_2,object['natoms_2']*g4b.perCent)
+        translation = g4b.G4ThreeVector(object['position_x']*g4b.um, object['position_y']*g4b.um, object['position_z']*g4b.um)
+        visual = g4b.G4VisAttributes(g4b.G4Color(object['colour'][0],object['colour'][1],object['colour'][2]))
+        mother = self.physical['world']
+        sidex = object['side_x']*g4b.um
+        sidey = object['side_y']*g4b.um
+        sidez = object['side_z']*g4b.um
+
+        if not 'tub' in object:
+            self.solid[name] = g4b.G4Box(name, sidex/2., sidey/2., sidez/2.)
+        else:
+            self.solid[name+"box"] = g4b.G4Box(name+"box", 
+                                           sidex/2., sidey/2., sidez/2.)
+            self.solid[name+"tub"] = g4b.G4Tubs(name+"tub", 0,object['tub']['tub_radius']*g4b.um,
+                                                object['tub']['tub_depth']*g4b.um, 0,360*g4b.deg)
+            self.solid[name] = g4b.G4SubtractionSolid(name,
+                                                    self.solid[name+"box"],
+                                                    self.solid[name+"tub"])
+            
+        self.logical[name] = g4b.G4LogicalVolume(self.solid[name], 
+                                                 compound, 
+                                                 name)
+        self.physical[name] = g4b.G4PVPlacement(None,translation,                                                
+                                                name,self.logical[name],
+                                                mother, False, 
+                                                0,self.checkOverlaps)
+        self.logical[name].SetVisAttributes(visual)
+
+    def Construct(self): 
+        return self.physical['world']
+
+class ActionInitialization(g4b.G4VUserActionInitialization):
+    def __init__(self, detConstruction, par_in, par_direction, par_type, par_energy, par_num):
+        super().__init__()
+        self.fDetConstruction = detConstruction
+        self.par_in = par_in
+        self.par_direction = par_direction
+        self.par_type = par_type
+        self.par_energy = par_energy
+        self.par_num = par_num
+
+        self.par_out = [0, 0, 0]
+
+    def Build(self):
+        self.SetUserAction(PrimaryGeneratorAction(self.par_in,
+                                                  self.par_direction,
+                                                  self.par_type,
+                                                  self.par_energy,
+                                                  self.par_num))
+                                                    
+        myRA_action = MyRunAction()
+        self.SetUserAction(myRA_action)
+        myEA = EventAction(myRA_action, self.par_in, self.par_out)
+        self.SetUserAction(myEA)
+        self.SetUserAction(MySteppingAction(myEA))
+
+class PrimaryGeneratorAction(g4b.G4VUserPrimaryGeneratorAction):
+    "My Primary Generator Action"
+    def __init__(self, par_in, par_direction, par_type, par_energy, par_num):
+        g4b.G4VUserPrimaryGeneratorAction.__init__(self)
+        self.par_num = par_num
+        self.par_direction = par_direction
+        self.par_energy = par_energy
+        self.par_in = par_in
+        particle_table = g4b.G4ParticleTable.GetParticleTable()
+        self.particle = particle_table.FindParticle(par_type) # define particle
+        
+    def GeneratePrimaries(self, event):
+        for i in range(self.par_num):
+            par_direction = self.par_direction[i]
+            par_in = self.par_in[i]
+            beam = g4b.G4ParticleGun(1)
+            beam.SetParticleEnergy(self.par_energy*g4b.MeV)
+            beam.SetParticleMomentumDirection(g4b.G4ThreeVector(par_direction[0],
+                                                                par_direction[1],
+                                                                par_direction[2]))
+            beam.SetParticleDefinition(self.particle)
+            beam.SetParticlePosition(g4b.G4ThreeVector(par_in[0]*g4b.um,
+                                                       par_in[1]*g4b.um,
+                                                       par_in[2]*g4b.um))
+            self.particleGun = beam
+
+            self.particleGun.GeneratePrimaryVertex(event)
+
+class MyRunAction(g4b.G4UserRunAction):
+    def __init__(self):
+        g4b.G4UserRunAction.__init__(self)
+      
+    def BeginOfRunAction(self, run):
+        g4b.G4RunManager.GetRunManager().SetRandomNumberStore(False)
+   
+    def EndOfRunAction(self, run):
+        nofEvents = run.GetNumberOfEvent()
+        if nofEvents == 0:
+            print("nofEvents=0")
+            return
+        
+class EventAction(g4b.G4UserEventAction):
+    "My Event Action"
+    def __init__(self, runAction, point_in, point_out):
+        g4b.G4UserEventAction.__init__(self)
+        self.fRunAction = runAction
+        self.point_in = point_in
+        self.point_out = point_out
+
+    def BeginOfEventAction(self, event):
+        self.edep_device=0.
+        self.event_angle = 0.
+        self.p_step = []
+        self.energy_step = []
+        #use in pixel_detector
+        self.volume_name = []
+        self.localposition = []
+        global eventID
+        eventID = event.GetEventID()
+
+    def EndOfEventAction(self, event):
+        eventID = event.GetEventID()
+        print("eventID:%s end"%eventID)
+        # if len(self.p_step):
+        #     point_a = [ b-a for a,b in zip(self.point_in,self.point_out)]
+        #     point_b = [ c-a for a,c in zip(self.point_in,self.p_step[-1])]
+        #     self.event_angle = cal_angle(point_a, point_b)
+        # else:
+        #     self.event_angle = None
+
+        # save_geant4_events(eventID, self.edep_device, self.p_step, self.energy_step, self.event_angle)
+        save_geant4_events(eventID, self.edep_device, self.p_step, self.energy_step)
+
+    def RecordDevice(self, edep, point_in, point_out):
+        self.edep_device += edep
+        self.p_step.append([point_in.getX()*1000,
+                           point_in.getY()*1000,point_in.getZ()*1000])
+        self.energy_step.append(edep)
+
+        
+def cal_angle(point_a,point_b):
+    "Calculate the angle between point a and b"
+    x=np.array(point_a)
+    y=np.array(point_b)
+    l_x=np.sqrt(x.dot(x))
+    l_y=np.sqrt(y.dot(y))
+    dot_product=x.dot(y)
+    if l_x*l_y > 0:
+        cos_angle_d=dot_product/(l_x*l_y)
+        angle_d=np.arccos(cos_angle_d)*180/np.pi
+    else:
+        angle_d=9999
+    return angle_d
+    
+def save_geant4_events(eventID, edep_device, p_step, energy_step):
+    print("save")
+    time.sleep(1)
+
+    if(len(p_step)>0):
+        s_eventIDs.append(eventID)
+        s_edep_devices.append(edep_device)
+        s_p_steps.append(p_step)
+        s_energy_steps.append(energy_step)
+        # s_events_angle.append(event_angle)
+    else:
+        s_eventIDs.append(eventID)
+        s_edep_devices.append(edep_device)
+        s_p_steps.append([[0,0,0]])
+        s_energy_steps.append([0])
+        # s_events_angle.append(event_angle)
+        
+class MySteppingAction(g4b.G4UserSteppingAction):
+    "My Stepping Action"
+    def __init__(self, eventAction):
+        g4b.G4UserSteppingAction.__init__(self)
+        self.fEventAction = eventAction
+
+    def UserSteppingAction(self, step):
+        edep = step.GetTotalEnergyDeposit()
+        point_pre  = step.GetPreStepPoint()
+        point_post = step.GetPostStepPoint() 
+        point_in   = point_pre.GetPosition()
+        point_out  = point_post.GetPosition()
+        volume = step.GetPreStepPoint().GetTouchable().GetVolume().GetLogicalVolume()
+        volume_name = volume.GetName()
+        # print(volume_name)
+        # time.sleep(0.5)
+
+        if(volume_name == "detector"): # important, no if => no signal
+            self.fEventAction.RecordDevice(edep, point_in, point_out)
\ No newline at end of file
diff --git a/bmos/get_signal.py b/bmos/get_signal.py
new file mode 100644
index 0000000000000000000000000000000000000000..e3b7d27df3186cab2b19833a77cdc2d8e39f33bd
--- /dev/null
+++ b/bmos/get_signal.py
@@ -0,0 +1,193 @@
+#!/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 get_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']
+
+#     print(my_d.device)
+#     print(voltage)
+    
+    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)
+
+    my_current = ccrt.CalCurrentG4P(my_d, my_f, my_g4p, 0)
+    totalengry = my_g4p.energy_steps
+
+    output_path = "raser/bmos/output/"
+    tag = f"{g4_dic['par_type']}_{g4_dic['par_energy']}MeV_{g4_dic['par_num']}particle"
+    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"
+
+    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)
+
+    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
+        # print(my_current.sum_cu[0][i])
+        # print(current_pwl)
+        # sleep(1)
+        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 read_file_current(file_path, file_name):
+    with open(os.path.join(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 draw(output_path, pwl_name, filename_after_ngspice, tag, totalengry):
+    file_path = output_path
+    
+    # geant4_json = "./setting/absorber/cflm.json"
+    # with open(geant4_json) as f:
+    #      g4_dic = json.load(f)
+
+    file_name_v = filename_after_ngspice
+    file_name_c = pwl_name
+
+    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(tag)
+    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(f"Energy Deposition:{sum(totalengry[0])}MeV")
+    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(os.path.join(output_path, f"Signal_{tag}.pdf"))
diff --git a/bmos/ucsc.cir b/bmos/ucsc.cir
new file mode 100644
index 0000000000000000000000000000000000000000..83dd11dea420d23d5709c1c1f79b8d97d9c3e6d4
--- /dev/null
+++ b/bmos/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/bmos/output/Current_test.raw v(out)
+* save raw data
+.endc
+
+.end
diff --git a/elec/set_pwl_input.py b/elec/set_pwl_input.py
index 9b4d0d4c3851949863a36cefc41263cf3bff761c..751d8e824dff7ecdfd6eb50ac0976e59bb6a3612 100644
--- a/elec/set_pwl_input.py
+++ b/elec/set_pwl_input.py
@@ -1,7 +1,7 @@
 import shutil
 import os
 
-def set_pwl_input(pwl_file, cir_file, output_folder):
+def set_pwl_input(pwl_file, cir_file, output_folder, output_filename = 'default'):
     
     
     string_list=[]
@@ -22,13 +22,20 @@ def set_pwl_input(pwl_file, cir_file, output_folder):
                 lines[i] = lines[i].strip().replace(' ', ',') + ','
                 string_list.append(lines[i])
     with open(new_cir_file, 'r') as f:
-         replacement_line = 'I1 2 0 PWL(' + ''.join(string_list) + ')'
-         output_lines = f.readlines()
-         for i in range(len(output_lines)):
-             if 'I1' in output_lines[i]:
-                 output_lines[i] = replacement_line + '\n'
+        replacement_line = 'I1 2 0 PWL(' + ''.join(string_list) + ')'
+        output_lines = f.readlines()
+        if output_filename == 'default':
+            for i in range(len(output_lines)):
+                if 'I1' in output_lines[i]:
+                    output_lines[i] = replacement_line + '\n'
+        else:
+            for i in range(len(output_lines)):
+                if 'I1' in output_lines[i]:
+                    output_lines[i] = replacement_line + '\n'
+                if 'wrdata' in output_lines[i]:
+                    output_lines[i] = f"wrdata {os.path.join(output_folder, output_filename)} v(out)" + '\n'
     with open(new_cir_file, 'w') as f:
-         f.writelines(output_lines)
+        f.writelines(output_lines)
     
     print('Temporary circuit file has been modified')
 
diff --git a/particle/carrier_list.py b/particle/carrier_list.py
index 120f4cf0d6e4a5c7f92553c6fc61f2486b0f22f4..d39b7bfbf271ca1730e2011881e3bac33a730ab9 100644
--- a/particle/carrier_list.py
+++ b/particle/carrier_list.py
@@ -20,7 +20,7 @@ class CarrierListFromG4P:
         elif (material == "Si"):
             self.energy_loss = 3.6 #ev
 
-        if batch == 0 and my_g4p.geant4_model == "time_resolution":
+        if batch == 0 and (my_g4p.geant4_model == "time_resolution" or my_g4p.geant4_model == "bmos"):
             total_step=0
             particle_number=0
             for p_step in my_g4p.p_steps_current:   # selecting particle with long enough track