Belle II Software development
generateRadBhabhas.py
1#!/usr/bin/env python
2
3
10
11# -------------------------------------------------------------------------------------
12# BG simulation: Radiative Bhabha's
13# usage:
14# basf2 generateRadBhabhas.py generator equivTime_us num [sampleType phase outdir]
15# arguments:
16# generator one of: bbbrem, bhwide, bhwide_largeangle
17# equivTime_us equivalent SuperKEKB running time in micro-seconds
18# num output file number
19# sampleType one of: study, usual, PXD, ECL (D = usual)
20# phase 2, 31 (= early phase 3, ie Run 1), 32 (= Run 2) or 3 (D = 3)
21# outdir output directory path (D = output)
22# -------------------------------------------------------------------------------------
23
24import basf2 as b2
25import sys
26import os
27from background import add_output
28
29# read parameters
30
31argvs = sys.argv
32argc = len(argvs)
33
34if argc == 4:
35 generator = argvs[1] # bbbrem, bhwide, bhwide_largeangle
36 equivTime = argvs[2] # equivalent SuperKEKB running time in micro-seconds
37 num = argvs[3] # output file number
38 sampleType = 'usual' # study, usual, PXD, ECL
39 phase = 3 # phase number
40 outdir = 'output' # output directory path
41elif argc == 5:
42 generator = argvs[1]
43 equivTime = argvs[2]
44 num = argvs[3]
45 sampleType = argvs[4]
46 phase = 3
47 outdir = 'output'
48elif argc == 6:
49 generator = argvs[1]
50 equivTime = argvs[2]
51 num = argvs[3]
52 sampleType = argvs[4]
53 phase = int(argvs[5])
54 outdir = 'output'
55elif argc == 7:
56 generator = argvs[1]
57 equivTime = argvs[2]
58 num = argvs[3]
59 sampleType = argvs[4]
60 phase = int(argvs[5])
61 outdir = argvs[6]
62else:
63 print('Usage:')
64 print(' basf2', argvs[0], 'generator equivTime_us num [sampleType phase outdir]')
65 print('Arguments:')
66 print(' generator one of: bbbrem, bhwide, bhwide_largeangle')
67 print(' equivTime_us equivalent SuperKEKB running time in micro-seconds')
68 print(' num output file number')
69 print(' sampleType one of: study, usual, PXD, ECL (D = usual)')
70 print(' phase 2, 31 (= early phase 3, ie Run 1), 32 (= Run 2) or 3 (D = 3)')
71 print(' outdir output directory path (D = output)')
72 sys.exit()
73
74# set parameters
75
76if generator == 'bbbrem':
77 bgType = 'RBB'
78 crossect = 524e6 * 2 # nb (factor of two: gamma emission by e- and by e+)
79elif generator == 'bhwide':
80 bgType = 'BHWide'
81 crossect = 123e3 # nb
82elif generator == 'bhwide_largeangle':
83 bgType = 'BHWideLargeAngle'
84 crossect = 123e3 # nb
85else:
86 print('unknown generator: ', generator)
87 sys.exit()
88
89if phase == 3:
90 lumi = 600 # /nb/s (1/nb/s = 1e33/cm^2/s)
91elif phase == 31:
92 lumi = 30 # /nb/s
93elif phase == 32:
94 lumi = 280 # /nb/s
95elif phase == 2:
96 lumi = 20 # /nb/s
97else:
98 print('phase ', phase, 'not supported')
99 sys.exit()
100
101realTime = float(equivTime) * 1000 # ns
102numEvents = int(crossect * lumi * realTime * 1e-9)
103fname = bgType + '_' + sampleType + '-phase' + str(phase) + '-' + num
104outputFile = outdir + '/' + fname + '.root'
105
106if numEvents == 0:
107 b2.B2ERROR('number of events is 0 -> increase equivTime_us')
108 sys.exit()
109
110# make output directory if it doesn't exist
111
112if not os.path.exists(outdir):
113 os.makedirs(outdir)
114
115# log messages
116
117b2.B2RESULT('Events to be generated: ' + str(numEvents) + ' - corresponds to ' + equivTime +
118 ' us of running at ' + str(lumi) + ' /nb/s (phase ' + str(phase) + ')')
119b2.B2RESULT('Output file: ' + outputFile)
120
121
122# Suppress messages and warnings during processing:
123b2.set_log_level(b2.LogLevel.RESULT)
124
125# Create path
126main = b2.create_path()
127kill = b2.create_path()
128
129
130def add_cut(name, minParticles, maxParticles, minTheta, maxTheta=None):
131 """Add a generator level cut and kill the event if the cut is not passed. In
132 this case the cut is on the min/max charged particles which have a
133 center-of-mass theta angle between minTheta and maxTheta. If maxTheta is not
134 given assume it to be 180-minTheta for a symmetric window"""
135
136 # if only one angle make it symmetric
137 if maxTheta is None:
138 maxTheta = 180 - minTheta
139 selection = main.add_module('GeneratorPreselection')
140 selection.param('applyInCMS', True)
141 selection.param('nChargedMin', minParticles)
142 selection.param('nChargedMax', maxParticles)
143 selection.param('MinChargedTheta', minTheta)
144 selection.param('MaxChargedTheta', maxTheta)
145 selection.param('MinChargedP', 0.)
146 selection.param('MinChargedPt', 0.)
147 selection.if_value("!=11", kill)
148 selection.set_name("generator cut: " + name)
149
150
151# Event info setter
152eventinfosetter = b2.register_module('EventInfoSetter')
153eventinfosetter.param('evtNumList', [numEvents])
154main.add_module(eventinfosetter)
155
156# Gearbox
157gearbox = b2.register_module('Gearbox')
158if phase == 2:
159 gearbox.param('fileName', 'geometry/Beast2_phase2.xml')
160elif phase == 31:
161 gearbox.param('fileName', 'geometry/Belle2_earlyPhase3.xml')
162elif phase == 32:
163 gearbox.param('fileName', 'geometry/Belle2_Run2.xml')
164if sampleType == 'study':
165 gearbox.param('override', [
166 ("/DetectorComponent[@name='PXD']//ActiveChips", 'true', ''),
167 ("/DetectorComponent[@name='PXD']//SeeNeutrons", 'true', ''),
168 ("/DetectorComponent[@name='SVD']//ActiveChips", 'true', ''),
169 ("/DetectorComponent[@name='SVD']//SeeNeutrons", 'true', ''),
170 ("/DetectorComponent[@name='TOP']//BeamBackgroundStudy", '1', ''),
171 ("/DetectorComponent[@name='ARICH']//BeamBackgroundStudy", '1', ''),
172 ("/DetectorComponent[@name='ECL']//BeamBackgroundStudy", '1', ''),
173 ("/DetectorComponent[@name='KLM']//BKLM/BeamBackgroundStudy", '1', ''),
174 ("/DetectorComponent[@name='KLM']//EKLM/BeamBackgroundStudy", '1', ''),
175 ])
176main.add_module(gearbox)
177
178# Generator
179if generator == "bbbrem":
180 main.add_module("BBBremInput", MinPhotonEnergyFraction=0.000001, Unweighted=True,
181 MaxWeight=1.57001e+07)
182 # at least one track below 0.5 degree means maximum one particle in 0.5-179.5
183 add_cut("at least one track below 0.5 degree", 0, 1, 0.5)
184elif generator == "bhwide":
185 main.add_module("BHWideInput", ScatteringAngleRangeElectron=[0.5, 179.5],
186 ScatteringAngleRangePositron=[0.5, 179.5])
187 add_cut("both tracks at least 0.5 degree", 2, 2, 0.5)
188 # but if one is above 1 and the other above 10 degree so we in 1-170 and
189 # 10-179
190 add_cut("max one track in 1-170", 0, 1, 1, 170)
191 add_cut("max one track in 10-179", 0, 1, 10, 179)
192elif generator == "bhwide_largeangle":
193 main.add_module("BHWideInput", ScatteringAngleRangeElectron=[0.5, 179.5],
194 ScatteringAngleRangePositron=[0.5, 179.5])
195 add_cut("both tracks at least 1 degree", 2, 2, 1)
196 add_cut("at least one 10 degree", 1, 2, 10)
197else:
198 print("unknown generation setting: {}".format(generator))
199
200# Geant geometry
201geometry = b2.register_module('Geometry')
202geometry.param('useDB', False)
203addComp = ["MagneticField3dQuadBeamline"]
204# add beast detectors
205if sampleType == 'study' and (phase == 31 or phase == 32):
206 addComp.extend(["BEAMABORT", "MICROTPC", "CLAWS", "HE3TUBE"])
207
208geometry.param({"excludedComponents": ["MagneticField"],
209 "additionalComponents": addComp})
210main.add_module(geometry)
211
212# Geant simulation
213fullsim = b2.register_module('FullSim')
214if sampleType == 'study':
215 fullsim.param('PhysicsList', 'FTFP_BERT_HP')
216 fullsim.param('UICommandsAtIdle', ['/process/inactivate nKiller'])
217 fullsim.param('StoreAllSecondaries', True)
218 fullsim.param('SecondariesEnergyCut', 0.000001) # [MeV] need for CDC EB neutron flux
219main.add_module(fullsim)
220
221# Show progress of processing
222progress = b2.register_module('Progress')
223main.add_module(progress)
224
225# Output
226if phase == 31 or phase == 32:
227 phase = 3
228add_output(main, bgType, realTime, sampleType, phase, fileName=outputFile)
229
230# Process events
231b2.process(main)
232
233# Print call statistics
234print(b2.statistics)