Belle II Software  release-08-01-10
RunSRBgMC.py
1 #!/usr/bin/env python3
2 
3 
10 
11 import basf2 as b2
12 import sys
13 argvs = sys.argv
14 argc = len(argvs)
15 
16 print("The input HEPEvt files needed to run SR MC sim are located in @kekcc:~igal/src/beast/SR")
17 
18 if len(sys.argv) != 4:
19  print("Usage: requires 2 arguments")
20  print("Argument 1: (SynchRad_HER | SynchRad_LER | test)")
21  print("Argument 2: file number")
22  print("Argument 3: ROOT output directory path")
23  sys.exit(1)
24 
25 name = argvs[1]
26 num = argvs[2]
27 output_dir = argvs[3]
28 # Set realTime you want to use
29 # realTime = 1.0e4 # 10us for each file <-- time in ns
30 # realTime = 20 # 20ns for each file <-- time in ns - TEST
31 
32 outputfilename = output_dir + "/output_" + name + "_" + num + ".root"
33 
34 # tagname = SynchRad_HER for HER and SynchRad_LER for LER
35 tagname = name
36 
37 
43 
44 
45 # suppress messages and warnings during processing:
46 b2.set_log_level(b2.LogLevel.ERROR)
47 # set_log_level(LogLevel.DEBUG)
48 
49 # PHASE 2
50 # FileIn = "Ph2_dt_4_8HER21445M.HEPEvt" # data for HER Phase2 dt_4-8 6.6um -> 1.07 of bunch current 0.8A
51 #
52 # FileIn = "Ph2_dt_4_8LER35124M.HEPEvt" #data for LERPhase2 dt_4-8
53 # 6.6um -> 1.4 of bunch current 1.0A
54 
55 eventinfosetter = b2.register_module('EventInfoSetter')
56 hepevtreader = b2.register_module('HepevtInput')
57 # Ph2_dt_4_8HER21445M.HEPEvt
58 # Ph2_dt_4_8HER21445MK2M.HEPEvt
59 # Ph2_dt_4_8LER35124M.HEPEvt
60 # PHASE 2 -> HER 6148 repetitions = 1ROF = 20us | LER 3560 repet = 1ROF = 20us
61 if name == "SynchRad_HER":
62  realTime = 1.0e4 # 10us per file
63  # data for HER Phase2 dt_4-8 6.6um -> 1.07 of bunch current 0.8A
64  FileIn = "Ph2_dt_4_8HER21445M.HEPEvt"
65  # 2340 - 10us realTime Ph2_dt_4_8HER21445M.HEPEvt
66  hepevtreader.param('inputFileList', [FileIn] * 2340)
67 elif name == "SynchRad_LER":
68  realTime = 1.0e4 # 10us per file
69  # data for LERPhase2 dt_4-8 6.6um -> 1.4 of bunch current 1.0A
70  FileIn = "Ph2_dt_4_8LER35124M.HEPEvt"
71  # 2340 - 10us realTime Ph2_dt_4_8HER21445M.HEPEvt
72  hepevtreader.param('inputFileList', [FileIn] * 2340)
73 elif name == "test":
74  realTime = 20 # 20ns per file
75  # data HER KeV -> MeV,is used for test to be sure that code works or for
76  # cross-check.
77  FileIn = "Ph2_dt_4_8HER21445MK2M.HEPEvt"
78  # - for quick TEST 5 ->~ 20nsec
79  hepevtreader.param('inputFileList', [FileIn] * 5)
80  name = "ynchRad_HER"
81 # hepevtreader.param('inputFileList', [FileIn]*1780) # 1780 - 10us realTime for Ph2_dt_4_8LER35124M.HEPEvt
82 # hepevtreader.param('inputFileList', [FileIn] * 2340) # 2340 - 10us
83 # realTime Ph2_dt_4_8HER21445M.HEPEvt
84 
85 # Register
86 gearbox = b2.register_module('Gearbox')
87 geometry = b2.register_module('Geometry')
88 simulation = b2.register_module('FullSim')
89 tagSetter = b2.register_module('BeamBkgTagSetter')
90 progress = b2.register_module('Progress')
91 
92 # Setting the option for all non-hepevt reader modules:
93 # number of events in the list must be >= number of entries in input file
94 # times number of repetitions
95 eventinfosetter.param('evtNumList', [1000000000])
96 eventinfosetter.param('runList', [1]) # from run number 1
97 eventinfosetter.param('expList', [1]) # and experiment number 1
98 
99 # PHASE 2
100 gearbox.param('fileName', '/geometry/Beast2_phase2.xml')
101 # select component you need, but always keep 'MagneticField', 'BeamPipe', 'Cryostat','HeavyMetalShield' !
102 # geometry.param('components', ['MagneticField', 'BeamPipe', 'Cryostat', 'PXD', 'HeavyMetalShield',
103 # 'SVD', 'BEAMABORT', 'PLUME', 'FANGS', 'CLAWS', 'CDC',
104 # 'ARICH', 'TOP', 'COIL', 'ECL', 'BKLM', 'EKLM', 'HE3', ''])
105 #
106 geometry.set_log_level(b2.LogLevel.INFO)
107 
108 # simulation.param('PhysicsList', "QGSP_BERT_HP")
109 simulation.param('PhysicsList', "QGSP_BERT_EMV") # faster than QGSP_BERT_HP
110 # simulation.param('PhysicsList', "FTFP_BERT_EMV")
111 simulation.param('UICommandsAtIdle', ["/process/inactivate nKiller"])
112 simulation.param("StoreAllSecondaries", True)
113 # in MeV we need this for CDC EB neutron flux simulation
114 simulation.param("SecondariesEnergyCut", 0.000001)
115 
116 main = b2.create_path()
117 
118 main.add_module(eventinfosetter)
119 main.add_module(progress)
120 main.add_module(hepevtreader)
121 # print_params(hepevtreader)
122 main.add_module(gearbox)
123 main.add_module(geometry)
124 main.add_module(simulation)
125 
126 # set background tag in SimHits, leave main path if all SimHits are empty
127 tagSetter.param('backgroundType', tagname)
128 tagSetter.param('realTime', realTime)
129 main.add_module(tagSetter)
130 emptyPath = b2.create_path()
131 tagSetter.if_false(emptyPath)
132 b2.print_params(tagSetter)
133 
134 # output: SimHits only
135 rootoutput = b2.register_module('RootOutput')
136 rootoutput.param('outputFileName', outputfilename)
137 rootoutput.param('updateFileCatalog', False)
138 
139 # PHASE 2
140 # Select branches you need in output.
141 rootoutput.param('branchNames',
142  ["SVDSimHits",
143  "SVDTrueHits",
144  "SVDTrueHitsToSVDSimHits",
145  "PXDSimHits",
146  "MCParticleToPXDSimHits",
147  "CLAWSSimHits",
148  "ClawsHits",
149  "FANGSSimHits",
150  "FANGSHits",
151  "PlumeSimHits",
152  "BeamabortSimHits",
153  "BeamabortHits",
154  "PindiodeSimHits",
155  "PindiodeHits",
156  "QcsmonitorSimHits",
157  "QcsmonitorHits",
158  "He3tubeSimHits",
159  "He3tubeHits",
160  "MicrotpcSimHits",
161  "MicrotpcHits",
162  "SADMetaHits"])
163 MIP_to_PE = [12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12]
164 he3digi = b2.register_module('He3Digitizer')
165 he3digi.param('conversionFactor', 0.303132019)
166 he3digi.param('useMCParticles', False)
167 main.add_module(he3digi)
168 diadigi = b2.register_module('BeamDigitizer')
169 diadigi.param('WorkFunction', 13.25)
170 diadigi.param('FanoFactor', 0.382)
171 main.add_module(diadigi)
172 pindigi = b2.register_module('PinDigitizer')
173 pindigi.param('WorkFunction', 3.64)
174 pindigi.param('FanoFactor', 0.13)
175 main.add_module(pindigi)
176 clawsdigi = b2.register_module('ClawsDigitizer')
177 clawsdigi.param('ScintCell', 16)
178 clawsdigi.param('C_keV_to_MIP', 457.114)
179 clawsdigi.param('C_MIP_to_PE', MIP_to_PE)
180 clawsdigi.param('PEthres', 1.0)
181 main.add_module(clawsdigi)
182 qcssdigi = b2.register_module('QcsmonitorDigitizer')
183 qcssdigi.param('ScintCell', 40)
184 qcssdigi.param('C_keV_to_MIP', 1629.827)
185 qcssdigi.param('C_MIP_to_PE', 15.0)
186 qcssdigi.param('MIPthres', 0.5)
187 main.add_module(qcssdigi)
188 fangsdigi = b2.register_module('FANGSDigitizer')
189 main.add_module(fangsdigi)
190 tpcdigi = b2.register_module('TpcDigitizer')
191 main.add_module(tpcdigi)
192 
193 main.add_module(rootoutput)
194 
195 # Process events
196 b2.process(main)
197 
198 # Print call statistics
199 print(b2.statistics)