fiedl
commented
2 years ago https://icecube-spno.slack.com/archives/D02ST977K6F/p1660137871514689
This should be a working example to propagate photons with current clsim:
#!/usr/bin/sh /cvmfs/icecube.opensciencegrid.org/py3-v4.2.1/icetray-start
#METAPROJECT /data/user/tyuan/temp/tarballs/icetray.main.r323f33be.Linux-x86_64.gcc-9.3.0
##METAPROJECT /data/user/tyuan/temp/tarballs/icetray.main.rf60f545a.Linux-x86_64.gcc-9.2.0
##METAPROJECT /scratch/tyuan/metaprojects/combo/build
import os
import argparse
import numpy as np
from I3Tray import I3Tray
from icecube import icetray, clsim, phys_services, simclasses, photonics_service
icetray.set_log_level(icetray.I3LogLevel.LOG_TRACE)
def main():
""" resim for monopod
"""
parser = argparse.ArgumentParser(
description='Resimulate benchmark i3s (e.g. under different ice)')
parser.add_argument('infiles', nargs='+')
parser.add_argument('-o', '--out', default='out.i3.zst',
help='output file')
parser.add_argument('-s', '--seed', default=8888, type=int,
help='seed for sprng')
parser.add_argument('-r', '--run', default=0, type=int,
help='run for sprng')
parser.add_argument('--nframes', type=int, default=None, help='number of frames to process')
parser.add_argument('--usecpus', default=False, action='store_true')
parser.add_argument('--domos', default=5, type=float)
parser.add_argument('--domeff', default=1., type=float)
parser.add_argument('--gcd', default=os.path.os.path.expandvars('$I3_DATA/GCD/GeoCalibDetectorStatus_2020.Run134142.Pass2_V0.i3.gz'),
type=str,
help='gcd file')
parser.add_argument('--icemodel', choices=('SpiceMie', 'SpiceLea', 'Spice3.2', 'Spice3.2.1',
'Spicebfr-v1', 'Spicebfr-v2', 'Spicebfr-v2_flat',
'Spicebfr-v2_tilt'),
default='Spicebfr-v2', help='The ice model to use for photon propagation.')
parser.add_argument('--usetables', default=False, action='store_true', help='use PhotonicsHitMaker')
parser.add_argument("--holeice",default='as.h2-50cm', choices=('as.flasher_p1_0.30_p2_0',
'as.flasher_p1_0.30_p2_-1',
'as.flasher_p1_0.30_p2_+1',
'as.flasher_p1_0.30_p2_-2',
'as.flasher_p1_0.30_p2_-3',
'as.h2-50cm',
'as.nominal',
'as.set0_p0=0.0_p1=0.0'), help="Holeice file")
parser.add_argument('--nocascadeextension', default=False, action='store_true')
parser.add_argument('--notilt', default=False, action='store_true')
parser.add_argument('--unweighted', default=False, action='store_true')
parser.add_argument('--scale', default=None, type=float)
parser.add_argument('--all', default=False, action='store_true')
parser.add_argument('--prescale', default=0.01, type=float)
parser.add_argument('--history', default=0, type=int)
args = parser.parse_args()
print("RNG being initiated...")
# set up a random number generator
randomService = phys_services.I3SPRNGRandomService(
seed = args.seed,
nstreams = 100000000,
streamnum = args.run)
tray = I3Tray()
tray.context['I3RandomService'] = randomService
tray.Add('I3Reader', Filenamelist=[args.gcd]+args.infiles)
tray.Add('Delete', keys=['InIcePulses',
'InIceDSTPulses',
'InIceRawData',
'IceTopDSTPulses',
'IceTopPulses',
'IceTopRawData',
'CleanIceTopRawData',
'ClusterCleaningExcludedTanks',
'OfflineIceTopHLCTankPulses',
'OfflineIceTopHLCVEMPulses',
'OfflineIceTopSLCVEMPulses',
'SimTrimmer',
'TankPulseMergerExcludedTanks',
'HLCPulsesTimeRange',
'I3SuperDST',
'TimeShift',
'BeaconLaunches',
'CalibratedSLC',
'CalibratedWaveformRange',
'CalibrationErrata',
'DSTTriggers',
'FilterMask_NullSplit0',
'PassedAnyFilter',
'PassedConventional',
'PassedKeepSuperDSTOnly',
'QFilterMask',
'I3MCPESeriesMap',
'I3MCPulseSeriesMap',
'I3MCPulseSeriesMapParticleIDMap',
'I3MCPESeriesMapWithoutNoise',
'PhotonSeriesMap',
'I3MCPESeriesMapParticleIDMap',
'SaturatedDOMs',
'SaturationWindows',
'UncleanedInIcePulsesTimeRange'])
if not args.usetables:
icebase = os.path.expandvars('$I3_SRC/ice-models/resources/models')
icemodel = os.path.join(icebase, 'ICEMODEL', f'{args.icemodel[:5].lower()}_{args.icemodel[5:].lower()}')
holeice = os.path.join(icebase, 'ANGSENS', f'angsens/{args.holeice}')
if not os.path.isfile(holeice):
holeice = os.path.join(icebase, 'ANGSENS', f'angsens_unified/{args.holeice}')
print(f'Running with {icemodel}, {holeice}')
# mostly from /data/user/jstachurska/HESE7/resimulation/photon_propagation/photon.py
# Note the automatic photon->MCPE conversion merges hits automatically as well :D
tray.Add(clsim.I3CLSimMakePhotons,
GCDFile=args.gcd,
UseCPUs=args.usecpus,
UseGPUs=not args.usecpus,
UseI3PropagatorService=False,
MCTreeName='I3MCTree',
OutputMCTreeName=None,
PhotonSeriesName='PhotonSeriesMap' if args.history > 0 or args.all else None,
MCPESeriesName=None if args.all else 'I3MCPESeriesMap',
RandomService=tray.context['I3RandomService'],
IceModelLocation=icemodel,
UseCascadeExtension=not args.nocascadeextension,
DisableTilt=args.notilt,
DoNotParallelize=False,
DOMOversizeFactor=args.domos,
DOMEfficiency=args.domeff,
HoleIceParameterization=holeice,
UnWeightedPhotons=args.unweighted,
UnWeightedPhotonsScalingFactor=args.scale,
StopDetectedPhotons=not args.all,
SaveAllPhotons=args.all,
SaveAllPhotonsPrescale=args.prescale,
PhotonHistoryEntries=args.history,
)
else:
spline_table_dir = os.path.expandvars('$I3_DATA/photon-tables/splines')
tilttabledir = os.path.expandvars('$I3_SRC/photonics-service/resources/tilt/')
eff_distance = os.path.expandvars(
f'$I3_DATA/photon-tables/splines/cascade_effectivedistance_spice_{args.icemodel[5:].lower()}_z20.eff.fits')
if args.icemodel == 'Spicebfr-v2':
cascade_service = photonics_service.I3PhotoSplineService(
amplitudetable = os.path.join(spline_table_dir,'cascade_single_spice_bfr-v2_flat_z20_a5.abs.fits'),
timingtable = os.path.join(spline_table_dir, 'cascade_single_spice_bfr-v2_flat_z20_a5.prob.fits'),
effectivedistancetable = eff_distance,
tiltTableDir=tilttabledir)
else:
raise RuntimeError(f'photons.py: Implement {args.icemodel} table-sim')
cascade_service.SetEfficiencies(args.domeff)
# split the MCTree into a cascade-only and a track-only version
tray.AddModule("I3MCTreeHybridSimulationSplitter")
tray.AddModule("I3PhotonicsHitMaker", "hitsFromTheTable",
CascadeService = cascade_service,
TrackService = None, # tracks are handled by clsim
Input = 'I3MCTreeCascades',
Output = 'I3MCPESeriesMap',
RandomService=tray.context['I3RandomService'],
)
## DEBUG
# tray.AddSegment(clsim.I3CLSimMakeHitsFromPhotons,"makePhotons",
# GCDFile=args.gcd,
# PhotonSeriesName="PhotonSeriesMap",
# MCPESeriesName="I3MCPESeriesMap",
# RandomService=tray.context['I3RandomService'],
# DOMOversizeFactor=args.domos,
# UnshadowedFraction=args.domeff,
# HoleIceParameterization=holeice,
# MergeHits=False,
# IceModelLocation=icemodel
# )
# tray.Add(lambda frame: print('nph:', len(frame['PhotonSeriesMap'])), Streams=[icetray.I3Frame.DAQ])
# def npe(frame):
# _ = 0
# for mod in frame['I3MCPESeriesMap'].keys():
# _ += np.sum([_.npe for _ in frame['I3MCPESeriesMap'][mod]])
# print('npe:', _)
# tray.AddModule(npe, Streams=[icetray.I3Frame.DAQ])
# Delete all MCPEs we're not operating on
def empty_mcpe(frame):
return len(frame['I3MCPESeriesMap'])>0
tray.AddModule(empty_mcpe, Streams=[icetray.I3Frame.DAQ],
If=lambda frame:frame.Has('I3MCPESeriesMap'))
tray.Add('I3Writer', filename=args.out,
streams=[icetray.I3Frame.DAQ])
if args.nframes is None:
tray.Execute()
else:
tray.Execute(args.nframes)
if __name__=='__main__':
main()
I'm trying to propagate some photons with clsim and visualize their path with steamshovel.
clsim does provide example scripts for that. However, they do fail when run against acurrent icetray build.
Maybe, the example scripts are just outdated.