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