3 from beamparameters
import add_beamparameters
9 from background
import add_output
11 d = datetime.datetime.today()
12 print(d.strftime(
'job start: %Y-%m-%d %H:%M:%S\n'))
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)")
27 generator = sys.argv[1].lower()
29 sampleType = sys.argv[3]
32 digitization = sys.argv[6]
37 outputfilename = output_dir +
'/output_phase_' + argvs[4] +
'_' + generator +
'_' + num +
'.root'
40 seed = str(1234567 + int(num))
43 kill = basf2.create_path()
44 main = basf2.create_path()
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"""
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)
69 beamparameters = add_beamparameters(main,
"Y4S")
70 beamparameters.param(
"smearVertex",
True)
71 beamparameters.param(
"generateCMS",
False)
73 if generator ==
"bbbrem":
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)
80 add_cut(
"at least one track below 0.5 degree", 0, 1, 0.5)
81 elif generator ==
"bhwide":
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)
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'
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":
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')
110 main.add_module(koralw)
111 elif generator ==
"aafh":
115 main.add_module(
"EventInfoSetter", expList=1, runList=1, evtNumList=evtnum)
116 aafh = register_module(
'AafhInput')
144 'suppressionLimits': [1e100] * 4,
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],
174 main.add_module(aafh)
176 print(
"unknown generation setting: {}".format(generator))
178 print(
'generator ', generator)
179 print(
'bgtype', bgType)
180 print(
'reaTime', realTime)
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',
''),
199 gearbox.param(
'override', [(
'/Global/length',
'40.0',
'm')])
201 gearbox.param(
'fileName',
'/geometry/Beast2_phase2.xml')
202 main.add_module(gearbox)
204 geometry = register_module(
'Geometry')
206 "excludedComponents": [
"MagneticField"],
207 "additionalComponents": [
"MagneticField3dQuadBeamline"],
210 main.add_module(geometry)
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)
222 main.add_module(fullsim)
224 main.add_module(
"Progress")
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",
235 "BeamabortSimHits",
"BeamabortHits",
236 "PindiodeSimHits",
"PindiodeHits",
237 "QcsmonitorSimHits",
"QcsmonitorHits",
238 "He3tubeSimHits",
"He3tubeHits",
239 "MicrotpcSimHits",
"MicrotpcHits",
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)
272 add_output(main, bgType, realTime, sampleType, phase, outputfilename)
277 print(
'Event Statistics:')
278 print(basf2.statistics)
280 d = datetime.datetime.today()
281 print(d.strftime(
'job finish: %Y-%m-%d %H:%M:%S\n'))