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