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