#!/usr/bin/env python # Copyright (c) 2019, Anthony Latorre # # 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 . """ 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)