diff options
| -rw-r--r-- | utils/mcpl | 78 |
1 files changed, 78 insertions, 0 deletions
diff --git a/utils/mcpl b/utils/mcpl new file mode 100644 index 0000000..5cfb949 --- /dev/null +++ b/utils/mcpl @@ -0,0 +1,78 @@ +#!/usr/bin/env python +# Copyright (c) 2019, Anthony Latorre <tlatorre at uchicago> +# +# This program is free software: you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the Free +# Software Foundation, either version 3 of the License, or (at your option) +# any later version. +# +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for +# more details. +# +# You should have received a copy of the GNU General Public License along with +# this program. If not, see <https://www.gnu.org/licenses/>. +""" +Script to recreate MCPL files for events which never got simulated in SNOMAN. +""" +from __future__ import print_function, division +import numpy as np +import pandas as pd +from sddm import read_hdf +import gzip + +def read_mcpl(filename): + data = np.genfromtxt(filename,skip_header=1) + rv = {} + n = 1 + i = 1 + with gzip.open(filename,"r") as f: + lines = f.readlines() + nevents = int(lines[0].split()[0]) + while True: + nparticles = int(data[i-1,0]) + rv[n] = lines[i:i+nparticles] + n += 1 + i += nparticles + + if i - 1 >= data.shape[0]: + break + + assert nevents == len(rv) + + return rv + +def get_mcgn(filename): + ev = read_hdf(filename, "ev") + mcgn = read_hdf(filename, "mcgn").assign(run=ev.run[0]) + mcgn = mcgn.groupby(['run','evn']).first().reset_index() + return mcgn + +if __name__ == '__main__': + import argparse + import sys + from os.path import join + + parser = argparse.ArgumentParser("recreate MCPL files") + parser.add_argument("filenames", nargs='+', help="atmospheric MC files") + parser.add_argument("--mcpl-directory", required=True, help="directory where MCPL files are stored") + parser.add_argument("-o", "--output", required=True, help="output directory") + args = parser.parse_args() + + for filename in args.filenames: + mcgn = get_mcgn(filename) + mcpl_filename = "mc_atm_genie_%06i_0.dat.gz" % mcgn.run[0] + mcpl_file = join(args.mcpl_directory,mcpl_filename) + mcpl = read_mcpl(mcpl_file) + for ev in mcpl.keys(): + if ev in mcgn.evn.values: + del mcpl[ev] + output = "mc_atm_genie_%06i_0.dat" % mcgn.run[0] + with open(join(args.output,output),"w") as fout: + fout.write("%i %i\n" % (len(mcpl), 19)) + for ev, lines in mcpl.iteritems(): + fout.write(''.join(lines)) + sys.exit(0) + + print(ev) |