import numpy as np
from NuRadioReco.utilities import units
from matplotlib import pyplot as plt
from radiotools import helper as hp
from NuRadioMC.utilities import cross_sections as cs
from NuRadioMC.utilities import earth_attenuation
from scipy.integrate import quad
from scipy.optimize import brentq
from scipy.interpolate import interp1d
from NuRadioMC.simulation.simulation import pretty_time_delta
from NuRadioMC.EvtGen.generator import write_events_to_hdf5
from NuRadioMC.utilities import inelasticities
import pickle
import os
import time
import logging
import NuRadioMC
from NuRadioReco.utilities import version
# np.random.seed(10)  # just for testing

logger = logging.getLogger(__name__)
logger.setLevel(logging.INFO)

R_earth = 6357390 * units.m
earth = earth_attenuation.PREM()


def generate_eventlist_cylinder(filename, n_events, Emin, Emax,
                                full_rmin=None, full_rmax=None, full_zmin=None, full_zmax=None,
                                thetamin=0.*units.rad, thetamax=np.pi * units.rad,
                                phimin=0.*units.rad, phimax=2 * np.pi * units.rad,
                                start_event_id=1,
                                flavor=[12, -12, 14, -14, 16, -16],
                                n_events_per_file=None,
                                spectrum='log_uniform',
                                start_file_id=0):
    """
    Event generator

    Generates neutrino interactions, i.e., vertex positions, neutrino directions,
    neutrino flavor, charged currend/neutral current and inelastiviy distributions.
    All events are saved in an hdf5 file.

    Parameters
    ----------
    filename: string
        the output filename of the hdf5 file
    n_events: int
        number of events to generate
    Emin: float
        the minimum neutrino energy (energies are randomly chosen assuming a
        uniform distribution in the logarithm of the energy)
    Emax: float
        the maximum neutrino energy (energies are randomly chosen assuming a
        uniform distribution in the logarithm of the energy)

    full_rmin: float (default None)
        lower r coordinate of simulated volume (if None it is set to 1/3 of the fiducial volume, if second vertices are not activated it is set to the fiducial volume)
    full_rmax: float (default None)
        upper r coordinate of simulated volume (if None it is set to 5x the fiducial volume, if second vertices are not activated it is set to the fiducial volume)
    full_zmin: float (default None)
        lower z coordinate of simulated volume (if None it is set to 1/3 of the fiducial volume, if second vertices are not activated it is set to the fiducial volume)
    full_zmax: float (default None)
        upper z coordinate of simulated volume (if None it is set to 5x the fiducial volume, if second vertices are not activated it is set to the fiducial volume)
    thetamin: float
        lower zenith angle for neutrino arrival direction
    thetamax: float
        upper zenith angle for neutrino arrival direction
    phimin: float
        lower azimuth angle for neutrino arrival direction
    phimax: float
         upper azimuth angle for neutrino arrival direction
    start_event: int
        default: 1
        event number of first event
    flavor: array of ints
        default: [12, -12, 14, -14, 16, -16]
        specify which neutrino flavors to generate. A uniform distribution of
        all specified flavors is assumed.
        The neutrino flavor (integer) encoded as using PDF numbering scheme,
        particles have positive sign, anti-particles have negative sign,
        relevant for us are:
        * 12: electron neutrino
        * 14: muon neutrino
        * 16: tau neutrino
    n_events_per_file: int or None
        the maximum number of events per output files. Default is None, which
        means that all events are saved in one file. If 'n_events_per_file' is
        smaller than 'n_events' the event list is split up into multiple files.
        This is useful to split up the computing on multiple cores.
    spectrum: string
        defines the probability distribution for which the neutrino energies are generated
        * 'log_uniform': uniformly distributed in the logarithm of energy
    start_file_id: int (default 0)
        in case the data set is distributed over several files, this number specifies the id of the first file
        (useful if an existing data set is extended)
        if True, generate deposited energies instead of primary neutrino energies
    """

    attributes = {}
    n_events = int(n_events)

    # save current NuRadioMC version as attribute
    # save NuRadioMC and NuRadioReco versions
    attributes['NuRadioMC_EvtGen_version'] = NuRadioMC.__version__
    attributes['NuRadioMC_EvtGen_version_hash'] = version.get_NuRadioMC_commit_hash()

    attributes['start_event_id'] = start_event_id

    attributes['fiducial_rmin'] = full_rmin
    attributes['fiducial_rmax'] = full_rmax
    attributes['fiducial_zmin'] = full_zmin
    attributes['fiducial_zmax'] = full_zmax
    attributes['rmin'] = full_rmin
    attributes['rmax'] = full_rmax
    attributes['zmin'] = full_zmin
    attributes['zmax'] = full_zmax
    attributes['flavors'] = flavor
    attributes['Emin'] = Emin
    attributes['Emax'] = Emax
    attributes['thetamin'] = thetamin
    attributes['thetamax'] = thetamax
    attributes['phimin'] = phimin
    attributes['phimax'] = phimax

    # define cylinder by two points and the radius
    h_cylinder = full_zmax - full_zmin
    r_cylinder = full_rmax
    pt1 = np.array([0, 0, R_earth + full_zmax])
    pt2 = np.array([0, 0, R_earth + full_zmin])

    data_sets = {}

    # calculate maximum width of projected area
    theta_max = np.arctan(h_cylinder / 2 / r_cylinder)
    d = 2 * r_cylinder * np.cos(theta_max) + h_cylinder * np.sin(theta_max)  # width of area

    print(f"cylinder r = {r_cylinder/units.km:.1f}km, h = {h_cylinder/units.km:.1f}km -> dmax = {d/units.km:.1f}km")

    def perp(a) :
        b = np.empty_like(a)
        b[0] = -a[1]
        b[1] = a[0]
        return b

    # line segment a given by endpoints a1, a2
    # line segment b given by endpoints b1, b2
    # return
    def seg_intersect(a1, a2, b1, b2) :
        da = a2 - a1
        db = b2 - b1
        dp = a1 - b1
        dap = perp(da)
        denom = np.dot(dap, db)
        num = np.dot(dap, dp)
        return (num / denom) * db + b1

    def get_R(t, v, X):
        """"
        calculate distance to center of Earth as a function of travel distance
        
        Parameters
        ----------
        t: 3dim array
            travel distance
        v: 3dim array
            direction
        X: 3dim array
            start point
        """
        return np.linalg.norm(v * t + X)

    def get_density(t, v, X):
        """
        calculates density as a function of travel distance
        
        Parameters
        ----------
        t: 3dim array
            travel distance
        v: 3dim array
            direction
        X: 3dim array
            start point
        """
        return  earth.density(get_R(t, v, X))

    def slant_depth(t, v, X):
        """
        calculates slant depth (grammage) as a function of travel distance
        
        Parameters
        ----------
        t: 3dim array
            travel distance
        v: 3dim array
            direction
        X: 3dim array
            start point
        """
        res = quad(get_density, 0, t, args=(v, X), limit=200)
        return res[0]

    def slant_depth_num(t, v, X, step=50 * units.m):
        """
        calculates slant depth (grammage) as a function of travel distance
        
        Parameters
        ----------
        t: 3dim array
            travel distance
        v: 3dim array
            direction
        X: 3dim array
            start point
        """
        tt = np.linspace(0, t, t / step)
        rr = np.linalg.norm(X + np.outer(tt, v), axis=1)
        res = np.trapz(earth.density(rr), tt)
        return res

    def obj_dist_to_surface(t, v, X):
        return get_R(t, v, X) - R_earth

    def obj(t, v, X, Lint):
        """
        objective function to determine at which travel distance we reached the interaction point
        """
        return slant_depth(t, v, X) - Lint

    def points_in_cylinder(pt1, pt2, r, q):
        """
        determines if point lies within a cylinder
        
        Parameters
        ----------
        pt1: 3dim array
            lowest point on cylinder axis
        pt2: 3dim array
            highest point on cylinder axis
        r: float
            radius of cylinder
        q: 3dim array
            point under test
        
        Returns True/False
        """
        vec = pt2 - pt1
        const = r * np.linalg.norm(vec)
        return len(np.where(np.dot(q - pt1, vec) >= 0 and np.dot(q - pt2, vec) <= 0 and np.linalg.norm(np.cross(q - pt1, vec)) <= const)[0]) > 0

    # precalculate the maximum slant depth to the detector
    if(not os.path.exists("buffer_Llimit.pkl")):
        zens = np.arange(0, 180.1 * units.deg, 2 * units.deg)
        Lint_max = np.zeros_like(zens)
        Lint_min = np.zeros_like(zens)
        Xs = np.array([[0, -0.5 * r_cylinder, -h_cylinder + R_earth],
                       [0, -0.5 * r_cylinder, R_earth],
                       [0, 0.5 * r_cylinder, -h_cylinder + R_earth],
                       [0, 0.5 * r_cylinder, R_earth]]
                       )
        v_tmps = -hp.spherical_to_cartesian(zens, np.zeros_like(zens))  # neutrino direction
        for i in range(len(v_tmps)):
            v = v_tmps[i]
            sdepth_tmp = np.zeros(4)
            for j, X in enumerate(Xs):
                if((X[2] == R_earth) and (zens[i] <= 90 * units.deg)):
                    t = 0
                    sdepth_tmp[j] = 0
                else:
                    t = brentq(obj_dist_to_surface, 100, 2 * R_earth, args=(-v, X))
                    sdepth_tmp[j] = slant_depth_num(t, -v, X)
        #         print(i, zens[i] / units.deg, X, sdepth_tmp[j])
        #     enter_point = X + (-v * t)
            Lint_max[i] = np.max(sdepth_tmp)
            Lint_min[i] = np.min(sdepth_tmp)
        pickle.dump([zens, Lint_max, Lint_min], open("buffer_Llimit.pkl", "wb"), protocol=4)
    else:
        zens, Lint_max, Lint_min = pickle.load(open("buffer_Llimit.pkl", "rb"))

    get_Lmax = interp1d(zens, Lint_max, kind='next')
    get_Lmin = interp1d(zens, Lint_min, kind='previous')

    if 0:
        fig, a = plt.subplots(1, 1)
        ztmp = np.linspace(0, 180 * units.deg, 10000)
        a.plot(ztmp / units.deg, get_Lmax(ztmp) / units.g * units.cm ** 2, 'C0-', label="max possible Lint")
    #     a.plot(zens / units.deg, Lint_max / units.g * units.cm ** 2, 'oC0')
        a.plot(ztmp / units.deg, get_Lmin(ztmp) / units.g * units.cm ** 2, 'C1-', label="min possible Lint")
    #     a.plot(zens / units.deg, Lint_min / units.g * units.cm ** 2, 'dC1')
        a.hlines(cs.get_interaction_length(.1 * units.EeV, 1, 12, "total") / units.g * units.cm ** 2, 0, 180, label="0.1 EeV", colors='C2')
        a.hlines(cs.get_interaction_length(1 * units.EeV, 1, 12, "total") / units.g * units.cm ** 2, 0, 180 , label="1 EeV", colors='C3')
        a.hlines(cs.get_interaction_length(10 * units.EeV, 1, 12, "total") / units.g * units.cm ** 2, 0, 180 , label="10 EeV", colors='C4')
        a.set_xlabel("zenith angle [deg]")
        a.set_ylabel("slant depth [g/cm^2]")
        a.semilogy(True)
        a.set_xticks(np.arange(0, 181, 10))
        a.set_ylim(5e5)
        a.legend()

        fig.tight_layout()
        fig.savefig("Lvszen.png")
        plt.show()

    failed = 0
    if(spectrum == 'log_uniform'):
        Enu = 10 ** np.random.uniform(np.log10(Emin), np.log10(Emax), n_events)
    flavors = np.array([flavor[i] for i in np.random.randint(0, high=len(flavor), size=n_events)])
    az = np.random.uniform(phimin, phimax, n_events)
    zen = np.arccos(np.random.uniform(np.cos(thetamax), np.cos(thetamin), n_events))
    # generate random positions on an area perpendicular do neutrino direction
    ax, ay = np.random.uniform(-0.5 * d, 0.5 * d, (2, n_events))
    # az = np.ones(n_events) * (R_earth - .5 * h_cylinder)  # move plane to the center of the cylinder

    # calculate grammage (g/cm^2) after which neutrino interacted
    Lint = np.random.exponential(cs.get_interaction_length(Enu, 1, flavors, "total"), n_events)

    mask = (Lint < get_Lmax(zen)) & (Lint > get_Lmin(zen))
    print(f"{np.sum(mask)}/{n_events} = {np.sum(mask)/n_events:.2g} can potentially interact in simulation volume")

    # calculate position where neutrino interacts
    data_sets = {'xx': [],
                'yy': [],
                'zz': [],
                'azimuths': [],
                'zeniths': [],
                'flavors': [],
                'energies': []}
    # calculate rotation matrix to transform position on area to 3D
    mask_int = np.zeros_like(mask, dtype=bool)
    t0 = time.perf_counter()
    n_cylinder = 0
    for j, i in enumerate(np.arange(n_events, dtype=int)[mask]):
        if(j % 1000 == 0 and i > 0):
            eta = (time.perf_counter() - t0) * (n_events - i) / i
            logger.info(f"{i}/{n_events} interacting = {np.sum(mask_int)}, failed = {failed}, n_cylinder = {n_cylinder}, eta = {pretty_time_delta(eta)}")
    #     print(f"calculating interaction point of event {i}"),
        c, s = np.cos(az[i]), np.sin(az[i])
        Raz = np.array(((c, -s, 0), (s, c, 0), (0, 0, 1)))
        c, s = np.cos(zen[i]), np.sin(zen[i])
    #     Rzen = np.array(((c, 1, -s), (0, 1, 0), (s, 0, c)))
        Rzen = hp.get_rotation(hp.spherical_to_cartesian(0, az[i]), hp.spherical_to_cartesian(zen[i], az[i]))
    #     R = hp.get_rotation(np.array([0, 0, 1]), hp.spherical_to_cartesian(zen[i], az[i]))
        R = np.matmul(Rzen, Raz)
        v = -hp.spherical_to_cartesian(zen[i], az[i])  # neutrino direction
        X = np.matmul(R, np.array([ax[i], ay[i], 0])) + np.array([0, 0, -0.5 * h_cylinder])
        if 0:
            import matplotlib.pyplot as plt
            from mpl_toolkits.mplot3d import Axes3D
            from mpl_toolkits.mplot3d.art3d import Poly3DCollection
            fig = plt.figure()
            a = fig.add_subplot(111, projection='3d')
            # Cylinder
            x = np.linspace(-r_cylinder, r_cylinder, 100)
            z = np.linspace(0, -h_cylinder, 100)
            Xc, Zc = np.meshgrid(x, z)
            Yc = np.sqrt(r_cylinder ** 2 - Xc ** 2)

            # Draw parameters
            rstride = 20
            cstride = 10
            a.plot_surface(Xc, Yc, Zc, alpha=0.2, rstride=rstride, cstride=cstride)
            a.plot_surface(Xc, -Yc, Zc, alpha=0.2, rstride=rstride, cstride=cstride)
            a.set_title(f"zenith = {zen[i]/units.deg:.0f}")
            xx = []
            yy = []
            zz = []
            for vert in np.array([[-0.5 * d, -0.5 * d, 0], [0.5 * d, -0.5 * d, 0], [0.5 * d, 0.5 * d, 0], [-0.5 * d, 0.5 * d, 0]]):
                t = np.matmul(Raz, vert) + np.array([0, 0, -0.5 * h_cylinder])
                xx.append(t[0])
                yy.append(t[1])
                zz.append(t[2])
            verts = [list(zip(xx, yy, zz))]
            a.add_collection3d(Poly3DCollection(verts, alpha=0.5))

            xx = []
            yy = []
            zz = []
            for vert in np.array([[-0.5 * d, -0.5 * d, 0], [0.5 * d, -0.5 * d, 0], [0.5 * d, 0.5 * d, 0], [-0.5 * d, 0.5 * d, 0]]):
                t = np.matmul(Rzen, np.matmul(Raz, vert)) + np.array([0, 0, -0.5 * h_cylinder])
                xx.append(t[0])
                yy.append(t[1])
                zz.append(t[2])
            verts = [list(zip(xx, yy, zz))]
            a.add_collection3d(Poly3DCollection(verts, alpha=0.5))

            s = np.array([0, 0, -0.5 * h_cylinder])
            t = v * 10 * units.km + s
            a.plot([s[0], t[0]], [s[1], t[1]], [s[2], t[2]], '-d')

            s = np.array([0, 0, -0.5 * h_cylinder])
            t = -hp.spherical_to_cartesian(90 * units.deg, az[i]) * 10 * units.km + s
            a.plot([s[0], t[0]], [s[1], t[1]], [s[2], t[2]], '--d')

            # check if neutrino axis is perpendicular
            t = np.array(verts[0][0]) - np.array(verts[0][1])
            t /= np.linalg.norm(t)
            print(np.dot(t, v))

            a.set_xlabel("x")
            a.set_zlabel("z")
            a.set_ylabel("y")
            plt.show()

        # check if trajectory passes through cylinder
    #     if(not points_in_cylinder(pt1, pt2, r_cylinder, X)):
        # we rotate everything in the plane defined by z and the propagration direction (such that v_y = 0)
        Xaz = np.matmul(Raz.T, X)
        rmin = Xaz[1]  # the closest distance to the z axis (center of cyllinder)
        if(abs(rmin) >= r_cylinder):
            continue
        # define the projected square of the cylinder
        # the two endpoints of the two horizontal lines are
        xtmp = (r_cylinder ** 2 - rmin ** 2) ** 0.5
        Lh1 = np.array([[-xtmp, 0], [xtmp, 0]])
        Lh2 = np.array([[-xtmp, -h_cylinder], [xtmp, -h_cylinder]])
        # the two endpoints of the two vertical lines are
        Lv1 = np.array([[-xtmp, 0], [-xtmp, -h_cylinder]])
        Lv2 = np.array([[xtmp, 0], [xtmp, -h_cylinder]])

        # define line of neutrino propagation by two points
        vaz = np.matmul(Raz.T, v)
        if(abs(vaz[1]) > 1e-10):
            a = 1 / 0
        v2d = np.array([vaz[0], vaz[2]])
        X2d = np.array([Xaz[0], Xaz[2]])
        t = 2 * d
        Paz = np.array([X2d + -t * v2d, X2d + t * v2d])

        # calculate points that intersect any of the 4 area (projected cylinder) boundaries
        intersects = []
        for k, (a1, a2) in enumerate([Lh1, Lh2]):
            tmp = seg_intersect(a1, a2, Paz[0], Paz[1])
            if((tmp[0] >= a1[0]) and (tmp[0] <= a2[0])):
                intersects.append(tmp)
        for k, (a1, a2) in enumerate([Lv1, Lv2]):
            tmp = seg_intersect(a1, a2, Paz[0], Paz[1])
            if((tmp[1] <= a1[1]) and (tmp[1] >= a2[1])):
                intersects.append(tmp)
        intersects = np.array(intersects)
        if(len(intersects) != 2):
            if 0:
                print(len(intersects))
                import matplotlib.pyplot as plt
                from mpl_toolkits.mplot3d import Axes3D
                fig = plt.figure()
                a = fig.add_subplot(111, projection='3d')
                # Cylinder
                x = np.linspace(-r_cylinder, r_cylinder, 100)
                z = np.linspace(0, -h_cylinder, 100)
                Xc, Zc = np.meshgrid(x, z)
                Yc = np.sqrt(r_cylinder ** 2 - Xc ** 2)

                # Draw parameters
                rstride = 20
                cstride = 10
                a.plot_surface(Xc, Yc, Zc, alpha=0.2, rstride=rstride, cstride=cstride)
                a.plot_surface(Xc, -Yc, Zc, alpha=0.2, rstride=rstride, cstride=cstride)
                X_enter = X + 10 * units.km * v
                X_leave = X - 10 * units.km * v
                a.plot([X_enter[0], X_leave[0]], [X_enter[1], X_leave[1]], [X_enter[2], X_leave[2]], '-o')

                a.set_xlabel("x")
                a.set_zlabel("z")
                a.set_ylabel("y")
                a.legend()
                a.set_title("no intersection")

                fig = plt.figure()
                a = fig.add_subplot(111, projection='3d')
                for (a1, a2) in [Lh1, Lh2, Lv1, Lv2]:
                    a.plot([a1[0], a2[0]], [rmin, rmin], [a1[1], a2[1]])
                a.plot([Paz[0][0], Paz[1][0]], [rmin, rmin], [Paz[0][1], Paz[1][1]])
                a.set_xlabel("x")
                a.set_zlabel("z")
                a.set_ylabel("y")
                a.legend()
                plt.show()
            continue  # neutrino is not passing through cylinder
        n_cylinder += 1
        ss = []
        for tmp in intersects:
            ss.append(np.dot(tmp - X2d, v2d.T))
        argsort = np.argsort(np.array(ss))  # check which intersection happens first along neutrino path
        if 0:
            if(len(intersects)):
                import matplotlib.pyplot as plt
                from mpl_toolkits.mplot3d import Axes3D
                fig = plt.figure()
                a = fig.add_subplot(111, projection='3d')
                for (a1, a2) in [Lh1, Lh2, Lv1, Lv2]:
                    a.plot([a1[0], a2[0]], [rmin, rmin], [a1[1], a2[1]])
                a.plot([Paz[0][0], Paz[1][0]], [rmin, rmin], [Paz[0][1], Paz[1][1]])
                for k, tmp in enumerate(intersects):
                    a.plot([tmp[0]], [rmin], [tmp[1]], 'o', label=f"s = {ss[k]:.0f}")
                a.set_xlabel("x")
                a.set_zlabel("z")
                a.set_ylabel("y")
                a.legend()
    #             plt.ion()
                plt.show()

        # calculate the 3D points where the neutrino enters/leaves the cylinder and transform to outside Earth
        X_enter = np.matmul(Raz, np.array([intersects[argsort][0][0], rmin, intersects[argsort][0][1]])) + np.array([0, 0, R_earth])
        X_leave = np.matmul(Raz, np.array([intersects[argsort][1][0], rmin, intersects[argsort][1][1]])) + np.array([0, 0, R_earth])
        X += np.array([0, 0, R_earth])

        # calculate point where neutrino enters Earth
        if(np.linalg.norm(X_enter) > R_earth):  # if enter point is outside of Earth (can happen because cylinder does not account for Earth curvature)
            if(np.linalg.norm(X_leave) > R_earth):  # check if leave point is also outside of Earth (can also happen because cylinder does not account for Earth curvature)
                continue
            t = brentq(obj_dist_to_surface, 0, 5 * d, args=(-v, X_leave))
            enter_point = X_leave + (-v * t)
            X_enter = enter_point  # define point where neutrino enters the cylinder as the point where it enters the Earth
        else:
            t = brentq(obj_dist_to_surface, 0, 2 * R_earth, args=(-v, X_enter))
            enter_point = X_enter + (-v * t)
    #     logger.debug(f"zen = {zen[i]/units.deg:.0f}deg, trajectory enters Earth at {enter_point[0]:.1f}, {enter_point[0]:.1f}, {enter_point[0]:.1f}. Dist to core = {np.linalg.norm(enter_point)/R_earth:.5f}, dist to (0,0,R) = {np.linalg.norm(enter_point - np.array([0,0,R_earth]))/R_earth:.4f}")

        # check if event interacts at all
        # calcualte slant depth to point of entering cylinder
        t = np.linalg.norm(enter_point - X_enter)
        slant_depth_min = slant_depth(t, v, enter_point)
        if(t == 0):
            slant_depth_min = 0
        # calculate slant depth through the cylinder
        s = np.linalg.norm(X_leave - X_enter)
        slant_depth_max = slant_depth(s, v, X_enter) + slant_depth_min  # full slant depth from outside Earth to point when it leaves the cylinder

        if 0:
            import matplotlib.pyplot as plt
            from mpl_toolkits.mplot3d import Axes3D
            fig = plt.figure()
            a = fig.add_subplot(111, projection='3d')
            # Cylinder
            x = np.linspace(-r_cylinder, r_cylinder, 100)
            z = np.linspace(0, -h_cylinder, 100)
            Xc, Zc = np.meshgrid(x, z)
            Yc = np.sqrt(r_cylinder ** 2 - Xc ** 2)

            # Draw parameters
            rstride = 20
            cstride = 10
            a.plot_surface(Xc, Yc, Zc, alpha=0.2, rstride=rstride, cstride=cstride)
            a.plot_surface(Xc, -Yc, Zc, alpha=0.2, rstride=rstride, cstride=cstride)
            a.plot([X_enter[0], X_leave[0]], [X_enter[1], X_leave[1]], [X_enter[2] - R_earth, X_leave[2] - R_earth], '-o')

            a.set_xlabel("x")
            a.set_zlabel("z")
            a.set_ylabel("y")
            print(f"Lmin = {slant_depth_min:.2g}, Lmax = {slant_depth_max:.2g}, Lnu = {Lint[i]:.2g}")
            a.legend()
            plt.show()
    #         a = 1 / 0
        if((Lint[i] <= slant_depth_min) or (Lint[i] >= slant_depth_max)):
            logger.debug("neutrino does not interact in cylinder, skipping to next event")
            continue

        try:
            # calculate interaction point by inegrating the density of Earth along the neutrino path until we wind the interaction length
            t = brentq(obj, 0, s, args=(v, X_enter, Lint[i] - slant_depth_min), maxiter=500)
        except:
            logger.warning("failed to converge, skipping event")
            failed += 1
            continue
        Xint = X_enter + v * t  # calculate interaction point

        is_in_cylinder = points_in_cylinder(pt1, pt2, r_cylinder, Xint)
        mask_int[i] = is_in_cylinder
        if(is_in_cylinder):
            logger.debug(f"event {i}, interaction point ({Xint[0]:.1f}, {Xint[1]:.1f}, {Xint[2]-R_earth:.1f}), in cylinder {is_in_cylinder}")
            data_sets['xx'].append(Xint[0])
            data_sets['yy'].append(Xint[1])
            data_sets['zz'].append(Xint[2] - R_earth)
            data_sets['zeniths'].append(zen[i])
            data_sets['azimuths'].append(az[i])
            data_sets['flavors'].append(flavors[i])
            data_sets['energies'].append(Enu[i])

        else:
            logger.error("interaction is not in cylinder but it should be")
            a = 1 / 0

    data_sets['event_ids'] = range(np.sum(mask_int))
    data_sets['flavors'] = np.ones(np.sum(mask_int))
    data_sets["event_ids"] = np.arange(np.sum(mask_int)) + start_event_id
    data_sets["n_interaction"] = np.ones(np.sum(mask_int), dtype=int)
    data_sets["vertex_times"] = np.zeros(np.sum(mask_int), dtype=float)

    data_sets["interaction_type"] = inelasticities.get_ccnc(np.sum(mask_int))
    data_sets["inelasticity"] = inelasticities.get_neutrino_inelasticity(np.sum(mask_int))

    attributes['n_events'] = n_events

    logger.info(f"{np.sum(mask_int)} event interacted in simulation volume")

    write_events_to_hdf5(filename, data_sets, attributes,
                         n_events_per_file=n_events_per_file, start_file_id=start_file_id)


if __name__ == "__main__":
    generate_eventlist_cylinder("test2.hdf5", 1e4, 1e18 * units.eV, 1e18 * units.eV, full_rmin=0, full_rmax=5 * units.km,
                                full_zmin=-2.7 * units.km, full_zmax=0, thetamin=0, thetamax=180 * units.deg,
                                phimin=0, phimax=360 * units.deg, start_event_id=0)
