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