Belle II Software  release-08-01-10
RunSADBgMC.py
1 #!/usr/bin/env python3
2 
3 
10 
11 import basf2 as b2
12 import sys
13 import datetime
14 from background import add_output
15 
16 d = datetime.datetime.today()
17 print(d.strftime('job start: %Y-%m-%d %H:%M:%S\n'))
18 
19 print("Make a link to the input SAD files before running this steering script:")
20 print("ln -s ~nakayama/[FIXME -> ask for the latest campaing] input")
21 print("OR")
22 print("Produce first the input SAD files before running this steering script:")
23 print("Make a link to the input SAD files before running this steering script:")
24 print("ln -s ~login_name/[FIXME -> ask your-self where did I put it] input")
25 
26 if len(sys.argv) != 9:
27  print("Usage: requires 8 arguments")
28  print("Argument 1: SAD input directory path")
29  print("Argument 2: (Touschek | Coulomb | Brens)_(HER | LER)(,_far)")
30  print("Argument 3: file number")
31  print("Argument 4: (study | usual | ECL | PXD)")
32  print("Argument 5: time-eqv in ns")
33  print("Argument 6: phase 1 | 2 |3")
34  print("Argument 7: ROOT output directory path")
35  print("Argument 8: digitization (true | false)")
36  sys.exit(1)
37 
38 # read parameters
39 argvs = sys.argv
40 argc = len(argvs)
41 
42 input_dir = argvs[1]
43 name = argvs[2]
44 num = argvs[3]
45 sampleType = argvs[4]
46 sampleTime = int(argvs[5])
47 phase = int(argvs[6])
48 output_dir = argvs[7]
49 digitization = argvs[8]
50 
51 print('Input dir: ', input_dir)
52 print('Name :', name)
53 print('File number: ', num)
54 print('Sample type: ', sampleType)
55 print('Sample time: ', sampleTime, 'ns')
56 print('Phase: ', phase)
57 print('Output dir', output_dir)
58 
59 # set accring (0:LER, 1:HER)
60 if name.find('LER') != -1:
61  accring = 0
62 elif name.find('HER') != -1:
63  accring = 1
64 else:
65  print('Name should include either of HER or LER')
66  sys.exit()
67 
68 if name.find('far') != -1:
69  range = 2800
70 else:
71  range = 400
72 
73 # inputfilename = input_dir + '/' + name + '.root'
74 fname = 'phase_' + str(argvs[6]) + '_' + name + '_' + sampleType + '_' + num
75 outputfilename = output_dir + '/output_' + fname + '.root'
76 bgType = name
77 
78 readouttime = 0
79 nevent = 1000000
80 # realTime = sampleTime # ns
81 
82 seed = str(1234567 + int(num))
83 nummod = str(int(num) % 5000)
84 
85 if phase == 2 or phase == 3:
86  realTime = 10.0e3 # ns
87  readmode = 0
88  inputfilename = 'input/EvtbyEvt/' + name + '_EvtbyEvt_' + nummod + '.root'
89  if name == 'Touschek_LER':
90  seed = seed + '3'
91  elif name == 'Touschek_HER':
92  seed = seed + '4'
93  elif name == 'Coulomb_LER':
94  seed = seed + '5'
95  elif name == 'Coulomb_HER':
96  seed = seed + '6'
97  elif name == 'Touschek_LER_far':
98  realTime = 0.1e3 # ns
99  readmode = 1
100  readouttime = 100 # 0.1us
101  seed = seed + '9'
102  elif name == 'Touschek_HER_far':
103  realTime = 0.1e3 # ns
104  readmode = 1
105  readouttime = 100 # 0.1us
106  seed = seed + '10'
107  elif name == 'Brems_HER':
108  seed = seed + '13'
109  elif name == 'Brems_LER':
110  readmode = 0
111  seed = seed + '14'
112  else:
113  print('Unknown name! (' + name + ')')
114  sys.exit()
115 elif phase == 1:
116  readmode = 1
117  readouttime = sampleTime
118  inputfilename = 'input/' + name + '_' + nummod + '.root'
119 else:
120  print('Unknown name! (' + name + ')')
121  sys.exit()
122 
123 if sampleType == 'study':
124  seed = seed + '1'
125 elif sampleType == 'usual':
126  seed = seed + '2'
127 elif sampleType == 'ECL':
128  seed = seed + '3'
129 elif sampleType == 'PXD':
130  seed = seed + '4'
131 else:
132  print('Unknown sample type! (' + sampleType + ')')
133  sys.exit()
134 
135 print('accring: ', accring, '(0:LER, 1:HER)')
136 print('input: ', inputfilename)
137 print('output: ', outputfilename)
138 print('range: ', range)
139 print('nevent: ', nevent)
140 print('readmode: ', readmode)
141 print('readouttime:', readouttime)
142 print('bgType: ', bgType)
143 print('sampleType: ', sampleType)
144 print('realTime: ', realTime, 'ns')
145 print('seed: ', seed)
146 
147 b2.set_log_level(b2.LogLevel.WARNING)
148 b2.set_random_seed(int(seed))
149 
150 main = b2.create_path()
151 
152 eventinfosetter = b2.register_module('EventInfoSetter')
153 eventinfosetter.param({'evtNumList': [nevent], 'runList': [1], 'expList': [1]})
154 main.add_module(eventinfosetter)
155 
156 gearbox = b2.register_module('Gearbox')
157 if sampleType == 'study' and phase == 3:
158  gearbox.param('override', [
159  ('/Global/length', '40.0', 'm'),
160  ("/DetectorComponent[@name='PXD']//ActiveChips", 'true', ''),
161  ("/DetectorComponent[@name='PXD']//SeeNeutrons", 'true', ''),
162  ("/DetectorComponent[@name='SVD']//ActiveChips", 'true', ''),
163  ("/DetectorComponent[@name='SVD']//SeeNeutrons", 'true', ''),
164  ("/DetectorComponent[@name='TOP']//BeamBackgroundStudy", '1', ''),
165  ("/DetectorComponent[@name='ARICH']//BeamBackgroundStudy", '1', ''),
166  ("/DetectorComponent[@name='ECL']//BeamBackgroundStudy", '1', ''),
167  ("/DetectorComponent[@name='BKLM']//BeamBackgroundStudy", '1', ''),
168  ("/DetectorComponent[@name='BeamPipe']//LimitStepLength", '1', ''),
169  ("/DetectorComponent[@name='Cryostat']//LimitStepLength", '1', ''),
170  ("/DetectorComponent[@name='FarBeamLine']//LimitStepLength", '1', ''),
171  ])
172 else:
173  gearbox.param('override', [('/Global/length', '40.0', 'm')]) # needed for FarBeamLine
174 if phase == 2:
175  gearbox.param('fileName', '/geometry/Beast2_phase2.xml')
176 elif phase == 1:
177  gearbox.param('fileName', '/geometry/Beast2_phase1.xml')
178 main.add_module(gearbox)
179 
180 geometry = b2.register_module('Geometry')
181 if phase == 2 or phase == 3:
182  geometry.param({
183  "excludedComponents": ["MagneticField"],
184  "additionalComponents": ["MagneticField3dQuadBeamline"],
185  })
186 main.add_module(geometry)
187 
188 sadinput = b2.register_module('SADInput')
189 sadinput.param('Filename', inputfilename)
190 sadinput.param('ReadMode', readmode)
191 sadinput.param('AccRing', accring)
192 sadinput.param('ReadoutTime', readouttime) # needed only for ReadMode = 1
193 sadinput.param('Range', range)
194 main.add_module(sadinput)
195 
196 fullsim = b2.register_module('FullSim')
197 fullsim.param('PhysicsList', 'FTFP_BERT_HP')
198 fullsim.param('UICommandsAtIdle', ['/process/inactivate nKiller'])
199 fullsim.param('StoreAllSecondaries', True)
200 fullsim.param('SecondariesEnergyCut', 0.0) # [MeV] need for CDC EB neutron flux
201 main.add_module(fullsim)
202 
203 progress = b2.register_module('Progress')
204 main.add_module(progress)
205 if phase == 1 and digitization == 'true':
206  rootoutput = b2.register_module('RootOutput')
207  rootoutput.param('outputFileName', outputfilename)
208  rootoutput.param('updateFileCatalog', False)
209  rootoutput.param('branchNames', ["ClawSimHits", "ClawsHits",
210  "BeamabortSimHits", "BeamabortHits",
211  "PindiodeSimHits", "PindiodeHits",
212  "BgoSimHits", "BgobortHits",
213  "CsiSimHits", "CsiHit_v2s",
214  "QcsmonitorSimHits", "QcsmonitorHits",
215  "He3tubeSimHits", "He3tubeHits",
216  "MicrotpcSimHits", "MicrotpcHits",
217  "SADMetaHits"])
218  bgodigi = b2.register_module('BgoDigitizer')
219  main.add_module(bgodigi)
220  # dosidigi = register_module('DosiDigitizer')
221  # main.add_module(dosidigi)
222  csidigi = b2.register_module('CsiDigitizer_v2')
223  main.add_module(csidigi)
224  he3digi = b2.register_module('He3Digitizer')
225  he3digi.param('conversionFactor', 0.303132019)
226  he3digi.param('useMCParticles', False)
227  main.add_module(he3digi)
228  diadigi = b2.register_module('BeamDigitizer')
229  diadigi.param('WorkFunction', 13.25)
230  diadigi.param('FanoFactor', 0.382)
231  main.add_module(diadigi)
232  pindigi = b2.register_module('PinDigitizer')
233  pindigi.param('WorkFunction', 3.64)
234  pindigi.param('FanoFactor', 0.13)
235  main.add_module(pindigi)
236  tpcdigi = b2.register_module('TpcDigitizer')
237  main.add_module(tpcdigi)
238  MIP_to_PE1 = [12.97, 12.46, 14.86, 15.71, 13.63, 14.56, 14.53, 15.31]
239  MIP_to_PE2 = [15.21, 12.46, 14.86, 15.71, 16.02, 15.83, 14.53, 15.31]
240  clawsdigi = b2.register_module('ClawDigitizer')
241  clawsdigi.param('ScintCell', 8)
242  clawsdigi.param('C_keV_to_MIP', 457.114)
243  clawsdigi.param('C_MIP_to_PE', MIP_to_PE2)
244  clawsdigi.param('PEthres', 1.0)
245  main.add_module(clawsdigi)
246  qcssdigi = b2.register_module('QcsmonitorDigitizer')
247  qcssdigi.param('ScintCell', 2)
248  qcssdigi.param('C_keV_to_MIP', 1629.827)
249  qcssdigi.param('C_MIP_to_PE', 15.0)
250  qcssdigi.param('MIPthres', 0.5)
251  main.add_module(qcssdigi)
252  main.add_module(rootoutput)
253 elif phase == 2 and digitization == 'true':
254  rootoutput = b2.register_module('RootOutput')
255  rootoutput.param('outputFileName', outputfilename)
256  rootoutput.param('updateFileCatalog', False)
257  rootoutput.param('branchNames', ["SVDSimHits", "SVDTrueHits", "SVDTrueHitsToSVDSimHits",
258  "PXDSimHits", "MCParticleToPXDSimHits",
259  "CLAWSSimHits", "ClawsHits",
260  "FANGSSimHits", "FANGSHits",
261  "PlumeSimHits",
262  "BeamabortSimHits", "BeamabortHits",
263  "PindiodeSimHits", "PindiodeHits",
264  "QcsmonitorSimHits", "QcsmonitorHits",
265  "He3tubeSimHits", "He3tubeHits",
266  "MicrotpcSimHits", "MicrotpcHits",
267  "SADMetaHits"])
268  MIP_to_PE = [12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12]
269  he3digi = b2.register_module('He3Digitizer')
270  he3digi.param('conversionFactor', 0.303132019)
271  he3digi.param('useMCParticles', False)
272  main.add_module(he3digi)
273  diadigi = b2.register_module('BeamDigitizer')
274  diadigi.param('WorkFunction', 13.25)
275  diadigi.param('FanoFactor', 0.382)
276  main.add_module(diadigi)
277  pindigi = b2.register_module('PinDigitizer')
278  pindigi.param('WorkFunction', 3.64)
279  pindigi.param('FanoFactor', 0.13)
280  main.add_module(pindigi)
281  clawsdigi = b2.register_module('ClawsDigitizer')
282  clawsdigi.param('ScintCell', 16)
283  clawsdigi.param('C_keV_to_MIP', 457.114)
284  clawsdigi.param('C_MIP_to_PE', MIP_to_PE)
285  clawsdigi.param('PEthres', 1.0)
286  main.add_module(clawsdigi)
287  qcssdigi = b2.register_module('QcsmonitorDigitizer')
288  qcssdigi.param('ScintCell', 40)
289  qcssdigi.param('C_keV_to_MIP', 1629.827)
290  qcssdigi.param('C_MIP_to_PE', 15.0)
291  qcssdigi.param('MIPthres', 0.5)
292  main.add_module(qcssdigi)
293  fangsdigi = b2.register_module('FANGSDigitizer')
294  main.add_module(fangsdigi)
295  tpcdigi = b2.register_module('TpcDigitizer')
296  main.add_module(tpcdigi)
297  main.add_module(rootoutput)
298 else:
299  add_output(main, bgType, realTime, sampleType, phase, outputfilename)
300 
301 b2.process(main)
302 
303 print('Event Statistics:')
304 print(b2.statistics)
305 
306 d = datetime.datetime.today()
307 print(d.strftime('job finish: %Y-%m-%d %H:%M:%S\n'))