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