Belle II Software  release-08-02-06
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 
24 import basf2 as b2
25 import sys
26 import os
27 from background import add_output
28 
29 # read parameters
30 
31 argvs = sys.argv
32 argc = len(argvs)
33 
34 if 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
41 elif argc == 5:
42  generator = argvs[1]
43  equivTime = argvs[2]
44  num = argvs[3]
45  sampleType = argvs[4]
46  phase = 3
47  outdir = 'output'
48 elif 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'
55 elif 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]
62 else:
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 
76 if generator == 'bbbrem':
77  bgType = 'RBB'
78  crossect = 524e6 * 2 # nb (factor of two: gamma emission by e- and by e+)
79 elif generator == 'bhwide':
80  bgType = 'BHWide'
81  crossect = 123e3 # nb
82 elif generator == 'bhwide_largeangle':
83  bgType = 'BHWideLargeAngle'
84  crossect = 123e3 # nb
85 else:
86  print('unknown generator: ', generator)
87  sys.exit()
88 
89 if phase == 3:
90  lumi = 600 # /nb/s (1/nb/s = 1e33/cm^2/s)
91 elif phase == 31:
92  lumi = 30 # /nb/s
93 elif phase == 32:
94  lumi = 280 # /nb/s
95 elif phase == 2:
96  lumi = 20 # /nb/s
97 else:
98  print('phase ', phase, 'not supported')
99  sys.exit()
100 
101 realTime = float(equivTime) * 1000 # ns
102 numEvents = int(crossect * lumi * realTime * 1e-9)
103 fname = bgType + '_' + sampleType + '-phase' + str(phase) + '-' + num
104 outputFile = outdir + '/' + fname + '.root'
105 
106 if 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 
112 if not os.path.exists(outdir):
113  os.makedirs(outdir)
114 
115 # log messages
116 
117 b2.B2RESULT('Events to be generated: ' + str(numEvents) + ' - corresponds to ' + equivTime +
118  ' us of running at ' + str(lumi) + ' /nb/s (phase ' + str(phase) + ')')
119 b2.B2RESULT('Output file: ' + outputFile)
120 
121 
122 # Suppress messages and warnings during processing:
123 b2.set_log_level(b2.LogLevel.RESULT)
124 
125 # Create path
126 main = b2.create_path()
127 kill = b2.create_path()
128 
129 
130 def 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
152 eventinfosetter = b2.register_module('EventInfoSetter')
153 eventinfosetter.param('evtNumList', [numEvents])
154 main.add_module(eventinfosetter)
155 
156 # Gearbox
157 gearbox = b2.register_module('Gearbox')
158 if phase == 2:
159  gearbox.param('fileName', 'geometry/Beast2_phase2.xml')
160 elif phase == 31:
161  gearbox.param('fileName', 'geometry/Belle2_earlyPhase3.xml')
162 elif phase == 32:
163  gearbox.param('fileName', 'geometry/Belle2_Run2.xml')
164 if 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  ])
176 main.add_module(gearbox)
177 
178 # Generator
179 if 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)
184 elif 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)
192 elif 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)
197 else:
198  print("unknown generation setting: {}".format(generator))
199 
200 # Geant geometry
201 geometry = b2.register_module('Geometry')
202 geometry.param('useDB', False)
203 addComp = ["MagneticField3dQuadBeamline"]
204 # add beast detectors
205 if sampleType == 'study' and (phase == 31 or phase == 32):
206  addComp.extend(["BEAMABORT", "MICROTPC", "CLAWS", "HE3TUBE"])
207 
208 geometry.param({"excludedComponents": ["MagneticField"],
209  "additionalComponents": addComp})
210 main.add_module(geometry)
211 
212 # Geant simulation
213 fullsim = b2.register_module('FullSim')
214 if 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
219 main.add_module(fullsim)
220 
221 # Show progress of processing
222 progress = b2.register_module('Progress')
223 main.add_module(progress)
224 
225 # Output
226 if phase == 31 or phase == 32:
227  phase = 3
228 add_output(main, bgType, realTime, sampleType, phase, fileName=outputFile)
229 
230 # Process events
231 b2.process(main)
232 
233 # Print call statistics
234 print(b2.statistics)