Belle II Software  release-08-01-10
RunLumiBgMC.py
1 
8 import basf2 as b2
9 from beamparameters import add_beamparameters
10 import sys
11 import datetime
12 from background import add_output
13 
14 d = datetime.datetime.today()
15 print(d.strftime('job start: %Y-%m-%d %H:%M:%S\n'))
16 
17 if len(sys.argv) != 7:
18  print("Usage: requires 5 arguments")
19  print("Argument 1: (bbbrem | bhwide | bhwide_largeangle | aafh | koralw)")
20  print("Argument 2: seed number")
21  print("Argument 3: (study | usual | ECL | PXD)")
22  print("Argument 4: phase ( 2 | 3)")
23  print("Argument 5: ROOT output directory path")
24  print("Argument 6: digitization (true | false)")
25  sys.exit(1)
26 
27 # read parameters
28 argvs = sys.argv
29 argc = len(argvs)
30 generator = sys.argv[1].lower()
31 num = argvs[2]
32 sampleType = sys.argv[3]
33 phase = int(argvs[4])
34 output_dir = argvs[5]
35 digitization = sys.argv[6]
36 
37 if phase == 2:
38  sampleType = 'usual'
39 
40 outputfilename = output_dir + '/output_phase_' + argvs[4] + '_' + generator + '_' + num + '.root'
41 
42 # set random seed
43 seed = str(1234567 + int(num))
44 b2.set_random_seed(seed)
45 
46 kill = b2.create_path()
47 main = b2.create_path()
48 # main.add_module("EventInfoSetter", expList=1, runList=1, evtNumList=1000000)
49 # main.add_module("EventInfoPrinter")
50 
51 # b2.set_log_level(b2.LogLevel.DEBUG)
52 # b2.set_debug_level(250)
53 # b2.logging.package("framework").log_level = b2.LogLevel.WARNING
54 
55 
56 def add_cut(name, minParticles, maxParticles, minTheta, maxTheta=None):
57  """Add a generator level cut and kill the event if the cut is not passed. In
58  this case the cut is on the min/max charged particles which have a
59  center-of-mass theta angle between minTheta and maxTheta. If maxTheta is not
60  given assume it to be 180-minTheta for a symmetric window"""
61 
62  # if only one angle make it symmetric
63  if maxTheta is None:
64  maxTheta = 180 - minTheta
65  selection = main.add_module("GeneratorPreselection", applyInCMS=True, nChargedMin=minParticles, nChargedMax=maxParticles,
66  MinChargedTheta=minTheta, MaxChargedTheta=maxTheta, MinChargedP=0., MinChargedPt=0.)
67  selection.if_value("!=11", kill)
68  selection.set_name("generator cut: " + name)
69 
70 
71 # beam parameters
72 beamparameters = add_beamparameters(main, "Y4S")
73 beamparameters.param("smearVertex", True)
74 beamparameters.param("generateCMS", False)
75 
76 if generator == "bbbrem":
77  bgType = 'RBB'
78  evtnum = 8430
79  realTime = 10 # 10ns
80  main.add_module("EventInfoSetter", expList=1, runList=1, evtNumList=evtnum)
81  main.add_module("BBBremInput", MinPhotonEnergyFraction=0.000001, Unweighted=True, MaxWeight=1.57001e+07)
82  # at least one track below 0.5 degree means maximum one particle in 0.5-179.5
83  add_cut("at least one track below 0.5 degree", 0, 1, 0.5)
84 elif generator == "bhwide":
85  bgType = 'BHWide'
86  evtnum = 984
87  realTime = 10.0e3 # 10us
88  main.add_module("EventInfoSetter", expList=1, runList=1, evtNumList=evtnum)
89  main.add_module("BHWideInput", ScatteringAngleRangeElectron=[0.5, 179.5], ScatteringAngleRangePositron=[0.5, 179.5])
90  add_cut("both tracks at least 0.5 degree", 2, 2, 0.5)
91  # but if one is above 1 and the other above 10 degree so we in 1-170 and
92  # 10-179
93  add_cut("max one track in 1-170", 0, 1, 1, 170)
94  add_cut("max one track in 10-179", 0, 1, 10, 179)
95 elif generator == "bhwide_largeangle":
96  bgType = 'BHWideLargeAngle'
97  evtnum = 98400
98  realTime = 1.0e6 # 1ms
99  main.add_module("EventInfoSetter", expList=1, runList=1, evtNumList=evtnum)
100  main.add_module("BHWideInput", ScatteringAngleRangeElectron=[0.5, 179.5], ScatteringAngleRangePositron=[0.5, 179.5])
101  add_cut("both tracks at least 1 degree", 2, 2, 1)
102  add_cut("at least one 10 degree", 1, 2, 10)
103 elif generator == "koralw":
104  bgType = 'twoPhoton'
105  evtnum = 1720
106  realTime = 1.0e3 # in ns = 1 us
107  main.add_module("EventInfoSetter", expList=1, runList=1, evtNumList=evtnum)
108  koralw = b2.register_module('KoralWInput')
109  koralw.param('RandomSeed', int(seed))
110  koralw.param('DataPath', './data/')
111  koralw.param('UserDataFile', 'KoralW_ee_mod.data')
112  # koralw.logging.log_level = LogLevel.INFO
113  main.add_module(koralw)
114 elif generator == "aafh":
115  bgType = 'twoPhoton'
116  evtnum = 5640
117  realTime = 1.0e3 # in ns = 1 us
118  main.add_module("EventInfoSetter", expList=1, runList=1, evtNumList=evtnum)
119  aafh = b2.register_module('AafhInput')
120  aafh.param({
121  # decay mode to generate.
122  # 1: e+e- -> mu+mu-L+L- where L is a user defined particle (default: tau)
123  # 2: e+e- -> mu+mu-mu+mu-
124  # 3: e+e- -> e+e-mu+mu-
125  # 4: e+e- -> e+e-L+L- where L is a user defined particle (default: tau)
126  # 5: e+e- -> e+e-e+e-
127  'mode': 5,
128  # to set the particle for modes 1 and 4 use set parameter "particle"
129  # rejection scheme to generate unweighted events
130  # 1: use rejection once for the final event weight
131  # 2: use rejection per sub generator and then for the final event
132  'rejection': 2,
133  # max subgenerator event weight, only used if rejection is set to 2
134  # (default). If this value is to low the generation will produce errors. If
135  # it is to high generation runs slower.
136  # 'maxSubgeneratorWeight': 1.0,
137  # max final event weight which is always used. If this value is to low the
138  # generation will produce errors. If it is to high generation runs slower.
139  # ==> should be around 2-4
140  # 'maxFinalWeight': 1.5,
141  # adjust subgenerator weights so that each sub generator has same
142  # probability to be called and the maximum weight is equal as well. These
143  # values are printed at the end of generation when output level is set to
144  # INFO. These weights strongly depend on the mode
145  # 'subgeneratorWeights': [1.000e+00, 2.216e+01, 3.301e+03, 6.606e+03, 1.000e+00, 1.675e+00, 5.948e+00, 6.513e+00],
146  # set to awfully precise
147  'suppressionLimits': [1e100] * 4,
148  # minimum invariant mass of the secondary pair
149  # 'minMass': 0.50,
150 
151  # 0.50GeV 2nd ###
152  # 'minMass': 0.50,
153  # 'maxFinalWeight': 3.59,
154  # 'maxSubgeneratorWeight': 0.716,
155  # 'subgeneratorWeights': [1.000e+00, 6.088e+01, 2.348e+04, 2.551e+04, 1.000e+00, 1.956e+00, 2.937e+00, 7.169e-01],
156 
157  # 0.01GeV 3rd ###
158  # 'minMass': 0.01,
159  # 'maxFinalWeight': 2.5,
160  # 'maxSubgeneratorWeight': 1.0,
161  # 'subgeneratorWeights': [1.000e+00, 2.216e+01, 3.301e+03, 6.606e+03, 1.000e+00, 1.675e+00, 5.948e+00, 6.513e+00],
162 
163  # 0.002GeV 1st ###
164  # 'minMass': 0.002,
165  # 'maxFinalWeight': 2.5,
166  # 'maxSubgeneratorWeight': 1.0,
167  # 'subgeneratorWeights': [1.000e+00, 2.216e+01, 3.301e+03, 6.606e+03, 1.000e+00, 1.675e+00, 5.948e+00, 6.513e+00],
168 
169  # 0.001GeV 1st ###
170  'minMass': 0.001,
171  'maxFinalWeight': 2.5,
172  'maxSubgeneratorWeight': 1.0,
173  'subgeneratorWeights': [1.000e+00, 2.216e+01, 3.301e+03, 6.606e+03, 1.000e+00, 1.675e+00, 5.948e+00, 6.513e+00],
174 
175  })
176  # aafh.logging.log_level = LogLevel.INFO
177  main.add_module(aafh)
178 else:
179  print("unknown generation setting: {}".format(generator))
180 
181 print('generator ', generator)
182 print('bgtype', bgType)
183 print('reaTime', realTime)
184 
185 gearbox = b2.register_module('Gearbox')
186 if sampleType == 'study' and phase == 3:
187  gearbox.param('override', [
188  ('/Global/length', '40.0', 'm'),
189  ("/DetectorComponent[@name='PXD']//ActiveChips", 'true', ''),
190  ("/DetectorComponent[@name='PXD']//SeeNeutrons", 'true', ''),
191  ("/DetectorComponent[@name='SVD']//ActiveChips", 'true', ''),
192  ("/DetectorComponent[@name='SVD']//SeeNeutrons", 'true', ''),
193  ("/DetectorComponent[@name='TOP']//BeamBackgroundStudy", '1', ''),
194  ("/DetectorComponent[@name='ARICH']//BeamBackgroundStudy", '1', ''),
195  ("/DetectorComponent[@name='ECL']//BeamBackgroundStudy", '1', ''),
196  ("/DetectorComponent[@name='BKLM']//BeamBackgroundStudy", '1', ''),
197  ("/DetectorComponent[@name='BeamPipe']//LimitStepLength", '1', ''),
198  ("/DetectorComponent[@name='Cryostat']//LimitStepLength", '1', ''),
199  ("/DetectorComponent[@name='FarBeamLine']//LimitStepLength", '1', ''),
200  ])
201 else:
202  gearbox.param('override', [('/Global/length', '40.0', 'm')]) # needed for FarBeamLine
203 if phase == 2:
204  gearbox.param('fileName', '/geometry/Beast2_phase2.xml')
205 main.add_module(gearbox)
206 
207 geometry = b2.register_module('Geometry')
208 geometry.param({
209  "excludedComponents": ["MagneticField"],
210  "additionalComponents": ["MagneticField3dQuadBeamline"],
211 })
212 # geometry.param('additionalComponents', ['FarBeamLine'])
213 main.add_module(geometry)
214 
215 # print generated particles
216 # mcparticleprinter = register_module('PrintMCParticles')
217 # mcparticleprinter.logging.log_level = LogLevel.INFO
218 # main.add_module(mcparticleprinter)
219 
220 fullsim = b2.register_module('FullSim')
221 fullsim.param('PhysicsList', 'FTFP_BERT_HP')
222 fullsim.param('UICommandsAtIdle', ['/process/inactivate nKiller'])
223 fullsim.param('StoreAllSecondaries', True)
224 fullsim.param('SecondariesEnergyCut', 0.0) # [MeV]
225 main.add_module(fullsim)
226 
227 main.add_module("Progress")
228 
229 if phase == 2 and digitization == 'true':
230  rootoutput = b2.register_module('RootOutput')
231  rootoutput.param('outputFileName', outputfilename)
232  rootoutput.param('updateFileCatalog', False)
233  rootoutput.param('branchNames', ["SVDSimHits", "SVDTrueHits", "SVDTrueHitsToSVDSimHits",
234  "PXDSimHits", "MCParticleToPXDSimHits",
235  "CLAWSSimHits", "ClawsHits",
236  "FANGSSimHits", "FANGSHits",
237  "PlumeSimHits",
238  "BeamabortSimHits", "BeamabortHits",
239  "PindiodeSimHits", "PindiodeHits",
240  "QcsmonitorSimHits", "QcsmonitorHits",
241  "He3tubeSimHits", "He3tubeHits",
242  "MicrotpcSimHits", "MicrotpcHits",
243  "SADMetaHits"])
244  MIP_to_PE = [12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12]
245  he3digi = b2.register_module('He3Digitizer')
246  he3digi.param('conversionFactor', 0.303132019)
247  he3digi.param('useMCParticles', False)
248  main.add_module(he3digi)
249  diadigi = b2.register_module('BeamDigitizer')
250  diadigi.param('WorkFunction', 13.25)
251  diadigi.param('FanoFactor', 0.382)
252  main.add_module(diadigi)
253  pindigi = b2.register_module('PinDigitizer')
254  pindigi.param('WorkFunction', 3.64)
255  pindigi.param('FanoFactor', 0.13)
256  main.add_module(pindigi)
257  clawsdigi = b2.register_module('ClawsDigitizer')
258  clawsdigi.param('ScintCell', 16)
259  clawsdigi.param('C_keV_to_MIP', 457.114)
260  clawsdigi.param('C_MIP_to_PE', MIP_to_PE)
261  clawsdigi.param('PEthres', 1.0)
262  main.add_module(clawsdigi)
263  qcssdigi = b2.register_module('QcsmonitorDigitizer')
264  qcssdigi.param('ScintCell', 40)
265  qcssdigi.param('C_keV_to_MIP', 1629.827)
266  qcssdigi.param('C_MIP_to_PE', 15.0)
267  qcssdigi.param('MIPthres', 0.5)
268  main.add_module(qcssdigi)
269  fangsdigi = b2.register_module('FANGSDigitizer')
270  main.add_module(fangsdigi)
271  tpcdigi = b2.register_module('TpcDigitizer')
272  main.add_module(tpcdigi)
273  main.add_module(rootoutput)
274 else:
275  add_output(main, bgType, realTime, sampleType, phase, outputfilename)
276 
277 # main.add_module("RootOutput", outputFileName="%s.root" % generator)
278 b2.process(main)
279 
280 print('Event Statistics:')
281 print(b2.statistics)
282 
283 d = datetime.datetime.today()
284 print(d.strftime('job finish: %Y-%m-%d %H:%M:%S\n'))