source/drich-dev/simulate.py

0001 #!/usr/bin/env python
0002
0003 # -----------------------------------------------#
0004 # npsim wrapper with EIC RICH specific tests     #
0005 # Author: C. Dilks                               #
0006 # -----------------------------------------------#
0007
0008 import sys, getopt, os, re, importlib
0009 import pprint
0010 import subprocess, shlex
0011 import math
0012 import numpy as np
0013
0014 ### convert theta <=> eta
0015 def theta_to_eta(th):
0016     return -math.log(math.tan(0.5 * th))
0017 def eta_to_theta(et):
0018     return 2 * math.atan(math.exp(-et))
0019
0020 # SETTINGS
0021 ################################################################
0022 acceptanceDict = {
0023     'drich': {
0024         # theta limits [degrees]
0025         'thetaMin': math.degrees(eta_to_theta(3.5)),
0026         'thetaMax': math.degrees(eta_to_theta(1.6)),
0027         # ideal theta: middle of acceptance, full rings in one sector
0028         'thetaIdeal': math.degrees(eta_to_theta(2.0)),
0029     },
0030     'pfrich': {
0031         # theta limits [degrees]
0032         'thetaMin': 180.0 - 10.0, # FIXME
0033         'thetaMax': 180.0 - 70.0, # FIXME
0034         # ideal theta: middle of acceptance
0035         'thetaIdeal': math.degrees(eta_to_theta(-2.0)), # FIXME
0036     },
0037 }
0038 momMax = {
0039   'aerogel': 20,
0040   'gas':     60,
0041 }
0042
0043 # ARGUMENTS
0044 ################################################################
0045
0046 ### defaults
0047 inputFileName = ''
0048 testNum = -1
0049 standalone = False
0050 compactFileCustom = ''
0051 zDirection = 1
0052 particle_name = 'pi+'
0053 particle_momentum = 20.0 # [GeV]
0054 particle_theta = acceptanceDict['drich']['thetaIdeal'] # FIXME: this default should be controllable by `zDirection`
0055 particle_eta = 10001
0056 runType = 'run'
0057 numEvents = 50
0058 numTestSamples = 0
0059 restrict_sector = True
0060 outputImageType = ''
0061 outputFileName = ''
0062
0063 ### usage guide
0064 helpStr = f'''
0065 {sys.argv[0]} <INPUT_FILE or TEST_NUM> [OPTIONS]
0066
0067 <REQUIRED ARGUMENTS>: provide either an INPUT_FILE or a TEST_NUM
0068
0069     INPUT_FILE: -i <input file>: specify an input file, e.g., hepmc
0070
0071     TEST_NUM:  -t <testnum>: specify which test to run
0072             >> acceptance tests:
0073
0074                fixed angle:
0075                 1:  aim pions at center of aerogel sector
0076                 2:  inner edge test
0077                 3:  outer edge test
0078
0079                acceptance scans: (-k = number of polar steps)
0080                                  (-n = number of particles per step)
0081                 4:  polar scan test, fixed azimuth
0082                 5:  azimuthal + polar scan test
0083                 6:  eta scan, varied azimuth for each eta value
0084
0085                momentum scans: (-k = number of momentum steps)
0086                                (-n = number of particles per step)
0087                 7:  momentum scan, for aerogel, fixed (theta,phi)
0088                 8:  momentum scan, for gas,     fixed (theta,phi)
0089                 9:  momentum scan, for aerogel, fixed theta, varied phi
0090                 10: momentum scan, for gas,     fixed theta, varied phi
0091
0092             >> optics tests: use these to test the dRICH optics;
0093                make sure to apply the recommended settings in drich.xml
0094                 100:   focal point, in RICH acceptance
0095                         set DRICH_debug_optics  = 1
0096                             DRICH_debug_mirror  = 0
0097                             DRICH_debug_sensors = 1
0098                 101:   focal point, broad range test
0099                         set DRICH_debug_optics  = 1
0100                             DRICH_debug_mirror  = 1
0101                             DRICH_debug_sensors = 1
0102                 102:   parallel-to-point focal test
0103                         set DRICH_debug_optics  = 1
0104                             DRICH_debug_mirror  = 0
0105                             DRICH_debug_sensors = 0
0106                 103:   evenly distributed sensor hits test
0107                         set DRICH_debug_optics  = 3
0108                             DRICH_debug_mirror  = 0
0109                             DRICH_debug_sensors = 0
0110                 104:   parallel-to-point focal test, beams over entire acceptance
0111                         set DRICH_debug_optics  = 4
0112                             DRICH_debug_mirror  = 0
0113                             DRICH_debug_sensors = 0
0114
0115 [OPTIONAL ARGUMENTS]
0116
0117     OPTIONS:    -d: direction to throw particles (may not be used by all tests)
0118                     1 = toward dRICH (default)
0119                    -1 = toward pfRICH
0120                 -s: enable standalone RICH-only simulation (default is full detector)
0121                 -c [compact file]: specify a custom compact file
0122                    (this will override -d and -s options)
0123                 -p [particle]: name of particle to throw; default: {particle_name}
0124                    examples:
0125                     - e- / e+
0126                     - pi+ / pi-
0127                     - kaon+ / kaon-
0128                     - proton / anti_proton
0129                     - opticalphoton
0130                 -m [momentum]: momentum (GeV) for mono-energetic runs (default={particle_momentum})
0131                 -a [angle]: fixed polar angle for certain tests [deg] (default={particle_theta})
0132                 -b [pseudorapidity]: fixed pseudorapidity for certain tests (default={theta_to_eta(math.radians(particle_theta))})
0133                      note: [pseudorapidity] will override [angle]
0134                 -n [numEvents]: number of events to process (default={numEvents})
0135                    - if using TEST_NUM, this is usually the number of events PER fixed momentum
0136                    - if using INPUT_FILE, you can set to 0 to run ALL events in the file, otherwise
0137                      it will run the default amount of {numEvents}
0138                 -k [numTestSamples]: some tests throw particles in multiple different directions,
0139                    such as "polar scan test"; for this test, use [numTestSamples] to control
0140                    how many directions are tested
0141                    - many tests offer a similar usage of [numTestSamples]
0142                    - these tests also have default [numTestSamples] values
0143                 -l: allow azimuthal scans to cover the full 2*pi range, rather than restricting
0144                     to a single sector
0145                 -r: run, instead of visualize (default)
0146                 -v: visualize, instead of run
0147                    - it is HIGHLY recommended to set `DRICH_debug_sector` to `1` in `drich.xml`,
0148                      which will draw one sector and set visibility such that you can see inside
0149                      the dRICH
0150                    - standalone mode will be automatically enabled
0151                 -e [output image extension]: save visual with specified type (svg,pdf,ps)
0152                    - useful tip: if you want to suppress the drawing of the visual, but
0153                      still save an output image, use Xvbf (start EIC container shell
0154                      as `xvfb-run eic-shell`); this is good for batch processing
0155                 -o [output file]: output root file name (overrides any default name)
0156     '''
0157
0158 if (len(sys.argv) <= 1):
0159     print(helpStr)
0160     sys.exit(2)
0161 try:
0162     opts, args = getopt.getopt(sys.argv[1:], 'i:t:d:sc:p:m:a:b:n:k:lrve:o:')
0163 except getopt.GetoptError:
0164     print('\n\nERROR: invalid argument\n', helpStr, file=sys.stderr)
0165     sys.exit(2)
0166 for opt, arg in opts:
0167     if (opt == '-i'): inputFileName = arg.lstrip()
0168     if (opt == '-t'): testNum = int(arg)
0169     if (opt == '-d'): zDirection = int(arg)
0170     if (opt == '-s'): standalone = True
0171     if (opt == '-c'): compactFileCustom = arg.lstrip()
0172     if (opt == '-p'): particle_name = arg.lstrip()
0173     if (opt == '-m'): particle_momentum = float(arg)
0174     if (opt == '-a'): particle_theta = float(arg)
0175     if (opt == '-b'): particle_eta = float(arg)
0176     if (opt == '-n'): numEvents = int(arg)
0177     if (opt == '-k'): numTestSamples = int(arg)
0178     if (opt == '-l'): restrict_sector = False
0179     if (opt == '-r'): runType = 'run'
0180     if (opt == '-v'): runType = 'vis'
0181     if (opt == '-e'): outputImageType = arg.lstrip()
0182     if (opt == '-o'): outputFileName = arg.lstrip()
0183 if (testNum < 0 and inputFileName == ''):
0184     print('\n\nERROR: Please specify either an input file (`-i`) or a test number (`-t`).\n', helpStr, file=sys.stderr)
0185     sys.exit(2)
0186 elif (testNum > 0 and inputFileName != ''):
0187     print('\n\nWARNING: You specified both an input file and a test number; proceeding with the input file only.\n', file=sys.stderr)
0188     testNum = -1
0189
0190 ### overrides
0191 if (testNum >= 100):
0192     print("optics test, overriding some settings...")
0193     particle_name = 'opticalphoton'
0194     standalone = True
0195     if (testNum in [100,101,102]):
0196         print("-- this is a visual test --")
0197         runType = 'vis'
0198 if (particle_name == "opticalphoton"):
0199     particle_momentum = 3e-9
0200     print(f'optical photons test: using energy {particle_momentum}')
0201 if runType == 'vis':
0202     standalone = True
0203
0204 ### set fixed particle angle & pseudorapidity: if particle_eta specified, override particle_theta
0205 if particle_eta < 10000:
0206     particle_theta = math.degrees(eta_to_theta(particle_eta))
0207 else:
0208     particle_eta = theta_to_eta(math.radians(particle_theta))
0209
0210 ### configure input and output file names
0211 ### relative paths will be made absolute here
0212 workDir = os.getcwd()
0213 ##### ensure input file name has absolute path
0214 if inputFileName != '':
0215     if not bool(re.search('^/', inputFileName)): inputFileName = workDir + "/" + inputFileName
0216 ##### ensure output file name has absolute path (and generate default name, if unspecified)
0217 if outputFileName == '':
0218     outputFileName = workDir + "/out/sim.edm4hep.root"  # default name
0219 elif not bool(re.search('^/', outputFileName)):
0220     outputFileName = workDir + "/" + outputFileName  # convert relative path to absolute path
0221 ##### get output file basename
0222 outputName = re.sub('\.root$', '', outputFileName)
0223 outputName = re.sub('^.*/', '', outputName)
0224
0225 ### set RICH names, based on zDirection
0226 zDirection /= abs(zDirection)
0227 if (zDirection < 0):
0228     xrich = 'pfrich'
0229     XRICH = 'PFRICH'
0230     xRICH = 'pfRICH'
0231 else:
0232     xrich = 'drich'
0233     XRICH = 'DRICH'
0234     xRICH = 'dRICH'
0235
0236 ### get env vars
0237
0238 detMain = os.environ['DETECTOR_CONFIG']
0239 detPath = os.environ['DETECTOR_PATH']
0240 outDir  = os.environ['DRICH_DEV'] + '/out'
0241
0242 ### set compact file
0243 compactFileFull = detPath + '/' + detMain + '.xml'
0244 compactFileRICH = detPath + '/' + detMain + '_' + xrich + '_only.xml'
0245 compactFile = compactFileRICH if standalone else compactFileFull
0246 if compactFileCustom != '':
0247     if not bool(re.search('^/', compactFileCustom)):
0248         compactFileCustom = workDir + "/" + compactFileCustom  # convert relative path to absolute path
0249     compactFile = compactFileCustom
0250
0251 ### print args and settings
0252 sep = '-' * 40
0253 print(f'''
0254 {sep}
0255 ** simulation args **
0256 inputFileName  = {inputFileName}
0257 testNum        = {testNum}
0258 particle       = {particle_name}
0259 particle_theta = {particle_theta} deg
0260 particle_eta   = {particle_eta}
0261 numEvents      = {numEvents}
0262 numTestSamples = {numTestSamples}
0263 runType        = {runType}
0264 direction      = toward {xRICH}
0265 outputFileName = {outputFileName}
0266 outputName     = {outputName}
0267 compactFile    = {compactFile}
0268 {sep}
0269 ''')
0270
0271 # SETTINGS AND CONFIGURATION
0272 ################################################################
0273
0274 ### start macro file
0275 m = open(workDir + "/macro/macro_" + outputName + ".mac", 'w+')
0276
0277 ### common settings
0278 m.write(f'/control/verbose 2\n')
0279 m.write(f'/run/initialize\n')
0280 # m.write(f'/run/useMaximumLogicalCores\n')
0281
0282 ### visual settings
0283 if (runType == 'vis'):
0284     m.write(f'/vis/open OGL 800x800-0+0\n')  # driver
0285     m.write(f'/vis/scene/create\n')
0286     m.write(f'/vis/scene/add/volume\n')
0287     m.write(f'/vis/scene/add/axes 0 0 0 1 m\n')
0288     m.write(f'/vis/scene/add/trajectories smooth\n')
0289     m.write(f'/vis/scene/add/hits\n')
0290     m.write(f'/vis/sceneHandler/attach\n')
0291     # m.write(f'/vis/viewer/set/viewpointThetaPhi 115 65\n') # angled view
0292     # m.write(f'/vis/viewer/set/viewpointThetaPhi 0 0\n') # front view
0293     m.write(f'/vis/viewer/set/viewpointThetaPhi -90 -89\n')  # top view
0294     # m.write(f'/vis/viewer/set/viewpointThetaPhi 90 0\n') # side view
0295     # m.write(f'/vis/viewer/zoom 0.5\n')
0296     m.write(f'/vis/viewer/set/style wireframe\n')
0297     m.write(f'/vis/modeling/trajectories/create/drawByCharge\n')
0298     m.write(f'/vis/modeling/trajectories/drawByCharge-0/setRGBA 0 0.8 0 0 1\n')
0299     m.write(f'/vis/modeling/trajectories/drawByCharge-0/setRGBA 1 0 0.5 0.5 1\n')
0300
0301 ### append particle info
0302 m.write(f'/gps/verbose 2\n')
0303 m.write(f'/gps/particle {particle_name}\n')
0304 m.write(f'/gps/number 1\n')
0305
0306 ### convert momentum to kinetic energy, for mono-energetic gun
0307 def momentum_to_kinetic_energy(p, part):
0308     # first get the mass
0309     mass = 0.0
0310     if bool(re.search('^e[+-]$',part)):
0311         mass = 0.000510999
0312     elif bool(re.search('^pi[+-]$',part)):
0313         mass = 0.13957
0314     elif bool(re.search('^kaon[+-]$',part)):
0315         mass = 0.493677
0316     elif bool(re.search('proton$',part)):
0317         mass = 0.938272
0318     elif (part == "opticalphoton"):
0319         mass = 0.0
0320     else:
0321         print(f'WARNING: mass for particle "{part}" needs to be added to simulate.py; assuming momentum==energy for now', file=sys.stderr)
0322     # then convert to energy
0323     en = math.sqrt( math.pow(p,2) + math.pow(mass,2) )
0324     kin_en = en - mass # total energy = kinetic energy + rest energy
0325     print(f'Momentum {p} GeV converted to Kinetic Energy {kin_en} GeV')
0326     return kin_en
0327 fixed_energy = momentum_to_kinetic_energy(particle_momentum, particle_name)
0328
0329 ### append source settings
0330 m.write(f'/gps/position 0 0 0 cm\n')
0331
0332 # ACCEPTANCE LIMITS
0333 ################################################################
0334
0335 ### set angular acceptance limits
0336 thetaMin = math.radians(acceptanceDict[xrich]['thetaMin'])
0337 thetaMax = math.radians(acceptanceDict[xrich]['thetaMax'])
0338 etaMin = theta_to_eta(thetaMax)
0339 etaMax = theta_to_eta(thetaMin)
0340 print(sep)
0341 print('** acceptance limits **')
0342 print(f'thetaMin = {math.degrees(thetaMin)} deg')
0343 print(f'thetaMax = {math.degrees(thetaMax)} deg')
0344 print(f'etaMin = {etaMin}')
0345 print(f'etaMax = {etaMax}')
0346 print(sep)
0347
0348 evnum = 0 # event number counter (for logging)
0349
0350 # TEST SETTINGS
0351 ######################################
0352
0353 ### `switch testNum:`
0354
0355 if testNum == 1:
0356     m.write(f'\n# aim at +x {xRICH} sector\n')
0357     x = math.sin(math.radians(particle_theta))
0358     y = 0.0
0359     z = math.cos(math.radians(particle_theta)) * zDirection
0360     m.write(f'/gps/direction {x} {y} {z}\n')
0361     m.write(f'/gps/ene/mono {fixed_energy} GeV\n')
0362     m.write(f'/run/beamOn {numEvents}\n')
0363
0364 elif testNum == 2:
0365     m.write(f'\n# inner edge of acceptance\n')
0366     x = math.sin(thetaMin)
0367     y = 0.0
0368     z = math.cos(thetaMin) * zDirection
0369     m.write(f'/gps/direction {x} {y} {z}\n')
0370     m.write(f'/gps/ene/mono {fixed_energy} GeV\n')
0371     m.write(f'/run/beamOn {numEvents}\n')
0372
0373 elif testNum == 3:
0374     m.write(f'\n# outer edge of acceptance\n')
0375     x = math.sin(thetaMax)
0376     y = 0.0
0377     z = math.cos(thetaMax) * zDirection
0378     m.write(f'/gps/direction {x} {y} {z}\n')
0379     m.write(f'/gps/ene/mono {fixed_energy} GeV\n')
0380     m.write(f'/run/beamOn {numEvents}\n')
0381
0382 elif testNum == 4:
0383     m.write(f'\n# polar scan test\n')
0384     numTheta = 4 if numTestSamples==0 else numTestSamples  # number of theta steps
0385     if (runType == "vis"):
0386         m.write(f'/vis/scene/endOfEventAction accumulate\n')
0387         m.write(f'/vis/scene/endOfRunAction accumulate\n')
0388     for theta in list(np.linspace(thetaMin, thetaMax, numTheta)):
0389         x = math.sin(theta)
0390         y = 0.0
0391         z = math.cos(theta) * zDirection
0392         m.write(f'/gps/direction {x} {y} {z}\n')
0393         m.write(f'/gps/ene/mono {fixed_energy} GeV\n')
0394         m.write(f'/run/beamOn {numEvents}\n')
0395         for _ in range(numEvents):
0396             print(f'evnum = {evnum}   theta = {math.degrees(theta)} deg   eta = {theta_to_eta(theta)}')
0397             evnum += 1
0398
0399 elif testNum == 5:
0400     m.write(f'\n# polar+azimuthal scan test\n')
0401     numTheta = 4 if numTestSamples==0 else numTestSamples # number of theta steps
0402     numPhi = 24  # number of phi steps, prefer even multiple of 6 (12,24,36) to check sector boundaries
0403     if (runType == "vis"):
0404         m.write(f'/vis/scene/endOfEventAction accumulate\n')
0405         m.write(f'/vis/scene/endOfRunAction accumulate\n')
0406     print(f'SET theta range to {math.degrees(thetaMin)} to {math.degrees(thetaMax)} deg')
0407     for theta in list(np.linspace(thetaMin, thetaMax, numTheta)):
0408         for phi in list(np.linspace(0, 2 * math.pi, numPhi, endpoint=False)):
0409             if restrict_sector and (math.pi / 6 < phi < (2 * math.pi - math.pi / 6)): continue  # restrict to one sector
0410             if (abs(phi) > 0.001 and abs(theta - thetaMin) < 0.001): continue  # allow only one ring at thetaMin
0411             x = math.sin(theta) * math.cos(phi)
0412             y = math.sin(theta) * math.sin(phi)
0413             z = math.cos(theta) * zDirection
0414             m.write(f'/gps/direction {x} {y} {z}\n')
0415             m.write(f'/gps/ene/mono {fixed_energy} GeV\n')
0416             m.write(f'/run/beamOn {numEvents}\n')
0417
0418 elif testNum == 6:
0419     m.write(f'\n# eta scan, with varied azimuth for each eta bin\n')
0420     numEtaPoints = 3 if numTestSamples==0 else numTestSamples # number of eta bins
0421     for eta in list(np.linspace(etaMin, etaMax, numEtaPoints)):
0422         m.write(f'\n# eta={eta}  theta={math.degrees(eta_to_theta(eta))} deg\n')
0423         for phi in list(2 * np.pi * np.random.random_sample(size=numEvents) - np.pi): # number of random phi values = `numEvents`
0424             x = math.sin(eta_to_theta(eta)) * math.cos(phi)
0425             y = math.sin(eta_to_theta(eta)) * math.sin(phi)
0426             z = math.cos(eta_to_theta(eta)) * zDirection
0427             m.write(f'/gps/direction {x} {y} {z}\n')
0428             m.write(f'/gps/ene/mono {fixed_energy} GeV\n')
0429             m.write(f'/run/beamOn 1\n')
0430
0431 elif testNum in [7,8]:
0432     rad = 'aerogel' if testNum==7 else 'gas'
0433     m.write(f'\n# momentum scan for {rad}, fixed theta, fixed phi\n')
0434     numMomPoints = 10 if numTestSamples==0 else numTestSamples # number of momenta
0435     x = math.sin(math.radians(particle_theta))
0436     y = 0.0
0437     z = math.cos(math.radians(particle_theta)) * zDirection
0438     m.write(f'/gps/direction {x} {y} {z}\n')
0439     for mom in list(np.linspace(1, momMax[rad], numMomPoints)):
0440         en = momentum_to_kinetic_energy(mom, particle_name)
0441         m.write(f'/gps/ene/mono {en} GeV\n')
0442         m.write(f'/run/beamOn {numEvents}\n')
0443
0444 elif testNum in [9,10]:
0445     rad = 'aerogel' if testNum==9 else 'gas'
0446     m.write(f'\n# momentum scan for {rad}, fixed theta, varied phi\n')
0447     numMomPoints = 10 if numTestSamples==0 else numTestSamples # number of momenta
0448     for mom in list(np.linspace(1, momMax[rad], numMomPoints)):
0449         en = momentum_to_kinetic_energy(mom, particle_name)
0450         for phi in list(2 * np.pi * np.random.random_sample(size=numEvents) - np.pi): # number of random phi values = `numEvents`
0451             x = math.sin(math.radians(particle_theta)) * math.cos(phi)
0452             y = math.sin(math.radians(particle_theta)) * math.sin(phi)
0453             z = math.cos(math.radians(particle_theta)) * zDirection
0454             m.write(f'/gps/direction {x} {y} {z}\n')
0455             m.write(f'/gps/ene/mono {en} GeV\n')
0456             m.write(f'/run/beamOn 1\n')
0457
0458 elif testNum == 100:
0459     m.write(f'\n# opticalphoton scan test, {xRICH} range\n')
0460     m.write(f'/vis/scene/endOfEventAction accumulate\n')
0461     m.write(f'/gps/pos/type Point\n')
0462     m.write(f'/gps/pos/radius 0.1 mm\n')
0463     m.write(f'/gps/ang/type iso\n')
0464     m.write(f'/gps/ang/mintheta {math.pi - thetaMax} rad\n')
0465     m.write(f'/gps/ang/maxtheta {math.pi - thetaMin} rad\n')
0466     m.write(f'/gps/ang/minphi {math.pi} rad\n')
0467     m.write(f'/gps/ang/maxphi {math.pi + 0.01} rad\n')
0468     m.write(f'/gps/ene/mono {fixed_energy} GeV\n')
0469     m.write(f'/run/beamOn {numEvents}\n')
0470
0471 elif testNum == 101:
0472     m.write(f'\n# opticalphoton scan test, broad range\n')
0473     m.write(f'/vis/scene/endOfEventAction accumulate\n')
0474     m.write(f'/gps/pos/type Point\n')
0475     m.write(f'/gps/pos/radius 0.1 mm\n')
0476     m.write(f'/gps/ang/type iso\n')
0477     m.write(f'/gps/ang/mintheta {math.pi / 2} rad\n')
0478     m.write(f'/gps/ang/maxtheta {math.pi - thetaMin} rad\n')
0479     m.write(f'/gps/ang/minphi {math.pi} rad\n')
0480     m.write(f'/gps/ang/maxphi {math.pi + 0.01} rad\n')
0481     m.write(f'/gps/ene/mono {fixed_energy} GeV\n')
0482     m.write(f'/run/beamOn {numEvents}\n')
0483
0484 elif testNum == 102:
0485     numTheta = 5 if numTestSamples==0 else numTestSamples # number of theta steps
0486     m.write(f'\n# opticalphoton parallel-to-point focusing\n')
0487     m.write(f'/vis/scene/endOfEventAction accumulate\n')
0488     m.write(f'/vis/scene/endOfRunAction accumulate\n')
0489     m.write(f'/gps/pos/type Beam\n')
0490     m.write(f'/gps/ang/type beam1d\n')
0491     for theta in list(np.linspace(thetaMin, thetaMax, numTheta)):
0492         x = math.sin(theta)
0493         y = 0.0
0494         z = math.cos(theta) * zDirection
0495         m.write(f'/gps/ang/rot1 -{z} {y} {x}\n') # different coordinate system...
0496         m.write(f'/gps/pos/rot1 -{z} {y} {x}\n')
0497         m.write(f'/gps/pos/halfx 16 cm\n')  # parallel beam width
0498         m.write(f'/gps/ene/mono {fixed_energy} GeV\n')
0499         m.write(f'/run/beamOn {numEvents}\n')
0500
0501 elif testNum == 103:
0502     m.write(f'\n# evenly distributed sensor hits test\n')
0503     if runType == "vis":
0504         m.write(f'/vis/scene/endOfEventAction accumulate\n')
0505         m.write(f'/vis/scene/endOfRunAction accumulate\n')
0506
0507     from scripts import createAngles
0508     num_rings = 120 if numTestSamples==0 else numTestSamples  # number of concentric rings, type=int
0509     hit_density = 80  # amount of photon hits for the smallest polar angle, type=int
0510     angles = createAngles.makeAngles(thetaMin, thetaMax, num_rings, hit_density)  # list of angles
0511
0512     print(f'SET theta range to {math.degrees(thetaMin)} to {math.degrees(thetaMax)} deg')
0513     for angle in angles:
0514         theta, phi = angle[0], angle[1]
0515         if restrict_sector and (math.pi / 6 < phi < (2 * math.pi - math.pi / 6)): continue  # restrict to one sector
0516         if abs(phi) > 0.001 and abs(theta - thetaMin) < 0.001: continue  # allow only one ring at thetaMin
0517         x = math.sin(theta) * math.cos(phi)
0518         y = math.sin(theta) * math.sin(phi)
0519         z = math.cos(theta) * zDirection
0520         m.write(f'/gps/direction {x} {y} {z}\n')
0521         m.write(f'/gps/ene/mono {fixed_energy} GeV\n')
0522         m.write(f'/run/beamOn {numEvents}\n')
0523
0524 elif testNum == 104:
0525     m.write(f'\n# opticalphoton parallel-to-point focusing, full coverage\n')
0526     #m.write(f'/vis/scene/endOfEventAction accumulate\n')
0527     #m.write(f'/vis/scene/endOfRunAction accumulate\n')
0528     m.write(f'/gps/pos/type Beam\n')
0529     m.write(f'/gps/ang/type beam1d\n')
0530     def makeBasicAngles(theta_min, theta_max, num_theta, num_phi):
0531         angles = []
0532         thetas = np.linspace(theta_min, theta_max, num=num_theta)
0533         phis = np.linspace(-math.pi/6, math.pi/6, num=num_phi)
0534         for i in range(num_phi):
0535             if phis[i] < 0:
0536                 phis[i] = phis[i]+2*np.pi
0537         for i in thetas:
0538             for j in phis:
0539                 angles.append(tuple((i,j)))
0540         return angles
0541     angles = makeBasicAngles(thetaMin, thetaMax, 50, 50)
0542
0543     for angle in angles:
0544         theta, phi = angle[0], angle[1]
0545         if math.pi / 6 < phi < (2 * math.pi - math.pi / 6): continue  # restrict to one sector
0546         x = math.sin(theta) * math.cos(phi)
0547         y = math.sin(theta) * math.sin(phi)
0548         z = math.cos(theta) * zDirection
0549         m.write(f'/gps/direction {x} {y} {z} \n')
0550         m.write(f'/gps/pos/halfx 16 cm\n')  # parallel beam width
0551         m.write(f'/gps/ene/mono {fixed_energy} GeV\n')
0552         m.write(f'/run/beamOn {numEvents}\n')
0553
0554 elif testNum > 0:
0555     print("ERROR: unknown test number\n", file=sys.stderr)
0556     m.close()
0557     sys.exit(2)
0558
0559 ### finalize
0560 if (runType == "vis"):
0561     m.write(f'/vis/viewer/flush\n')
0562     m.write(f'/vis/viewer/refresh\n')
0563     if outputImageType!='':
0564         m.write(f'/vis/ogl/export {re.sub("root$",outputImageType,outputFileName)}\n')
0565
0566 ### print macro and close stream
0567 m.seek(0, 0)
0568 if (testNum > 0):
0569     print(m.read())
0570 m.close()
0571
0572 # RUN npsim
0573 #########################################################
0574
0575 ### simulation executable and arguments
0576 cmd = [
0577         f'npsim',
0578         # f'{localDir}/NPDet/install/bin/npsim', # call local npsim
0579         f'--runType {runType}',
0580         f'--compactFile {compactFile}',
0581         f'--outputFile {outputFileName}',
0582         "--part.userParticleHandler=''", # necessary for opticalphotons truth output
0583         # '--random.seed 1',
0584         # '--part.keepAllParticles True',
0585         ]
0586 if (testNum > 0):
0587     cmd.extend([
0588         f'--macro {m.name}',
0589         '--enableG4GPS',
0590         ])
0591 else:
0592     cmd.extend([
0593         f'--inputFiles \'{inputFileName}\'',
0594         ])
0595     if (numEvents > 0):
0596         cmd.extend([ f'-N {numEvents}' ])
0597     else:
0598         cmd.extend([ f'-N -1' ])
0599
0600 ### run simulation
0601 cmdShell = shlex.split(" ".join(cmd))
0602 print(f'{sep}\nRUN SIMULATION:\n{shlex.join(cmdShell)}\n{sep}')
0603 subprocess.run(cmdShell, cwd=detPath)
0604
0605 ### cleanup
0606 # os.remove(m.name) # remove macro
0607 print("\nPRODUCED SIMULATION OUTPUT FILE: " + outputFileName)