<jupyter_start><jupyter_text>Splotter **First, make sure you have these packages:**- numpy- pandas- matplotlib- astropy- linetools **Then, run the cell below and go to the last cell in this notebook under 'Example Code' and adjust the variables.****Make sure you have the files I sent on Slack downloaded to the directory you are working in.**<jupyter_code>def splotter(high, low1, low2, user_dir, spectra_dir, all_csv, galaxy_csv, diff_tol=40):
    '''
    Inputs:
    high: string containing name of high-ion (needs to match transition name in all_df dataframe)
    low1: string containing name of first low-ion (needs to match transition name in all_df dataframe)
    low2: string containing name of second low-ion (needs to match transition name in all_df dataframe)
    all_csv: csv containing including any high-ion and low-ions to be aligned in stack plots
    galaxy_csv: csv containing all galaxies. Any ions specified as high, low1 or low2 within 500km/s will be included in stack plots
    user_dir: string containing directory in which the user wish to save stack plots
    spectra_dir: string containing directory where the analyzed spectra (.fits files) for each sightline is located
    diff_tol: float that determines within what central velocity difference a line is not considered no-low (in km/s) (default is 40)
    
    Outputs:
    Saves pdfs (stack plots) for all alignments within the high-ion and low-ion dataframe
    These pdfs are saved in a directory named identically after the high-ion input string and in organized subdirectories categorized as "stack_plots", "stack_plots_data" and "stats"
    Both "stack_plots" and "stack_plots_data" have subdirectories "no-lows", "narrows" and "broads" corresponding to the identified alignment categories and all a copy of these plots are also saved to the "all_alignments" subdirectory
    Each plot/pdf that is saved to "stack_plots" has a corresponding .csv file which contains important data about each ion and the nearby galaxy and are saved to the subdirectory "stack_plots_data"
    
    '''
    
    # Importing necessary libraries:
    
    import os
    
    import numpy as np
    import pandas as pd
    import math
    import glob

    import matplotlib.pyplot as plt

    from astropy.io import fits
    from astropy.table import Table

    from astropy import units as u
    from astropy import constants as const
    # New definition
    c_in_km_per_s = const.c.to('km/s').value
    
    from matplotlib.ticker import (MultipleLocator, AutoMinorLocator)
    
    from linetools.spectralline import AbsLine
    from linetools.isgm import utils as ltiu
    from linetools.analysis import absline as laa
    from linetools.spectra import io as lsio
    from linetools.isgm.abscomponent import AbsComponent
    from linetools.spectra.xspectrum1d import XSpectrum1D

    import imp
    lt_path = imp.find_module('linetools')[1]

    from linetools.analysis import voigt as lav

    #new imports
    from linetools.lists.linelist import LineList
    from linetools.analysis.plots import stack_plot


    ISM_linelist = LineList('ISM')
    plt.style.use('fast')
    
    
    # Defining preliminary functions:
    
    def make_absline(df_row):
        restwave = df_row.restwave*u.AA
        z_absorber = (1 + df_row.zsys)*(1 + df_row.vel/c_in_km_per_s) - 1
        attribs = {
            'b': df_row.bval * u.km/u.s,
            'logN': df_row.col,
            'QSO': df_row.QSO,
            'N':10**df_row.col / u.cm**2,
            'vel': df_row.vel * u.km/u.s,
            'z_abs': z_absorber,
            'z_gal': df_row.z_gal,
            'trans': df_row.trans
                  } 
        #lav.voigt_from_abslines wants the column density info to be in the `'N'` attrib, doesn't check `'logN'`

        #pass in a linelist so `linetools` doesn't load one every time this method is called
        absline = AbsLine(restwave, z=z_absorber, linelist=ISM_linelist)
        for name, value in attribs.items():
            absline.attrib[name] = value

        return absline
    
    def splotter_plotter(high, low1, low2, lines, aligned_df, galaxy_df, qso, z_high, keyword, user_dir, spectra_dir): #subdir takes an array of strings
        if keyword == "no-low":
            subdirs = ['/no-lows', '/all_alignments']
        elif keyword == "narrow":
            subdirs = ['/narrows', '/all_alignments']
        elif keyword == "broad":
            subdirs = ['/broads', '/all_alignments']
        else:
            print('INVALID KEYWORD')


        for index in range(0, len(subdirs)):
            #`indices` are all associated with the same QSO and the system
            z_galaxy = z_high
            sightline = qso #line.attrib['QSO']
            xspec = lsio.readspec(f'{spectra_dir}{sightline}_nbin3_norm.fits')

            for line in lines:
                line.analy['spec'] = xspec # Start plotting

            #synthesize and assign model voigt profiles
            wv_array = xspec.wavelength
            all_voigts = lav.voigt_from_abslines(wv_array, lines, fwhm=None, 
                                         ret=['flux'], debug=False)

            #all_voigts.sig = np.full([len(wv_array)], 0)
            xspec.sig = all_voigts * xspec.co

            fig = stack_plot(lines, zref=z_galaxy, vlim=[-750, 750]*u.km/u.s, return_fig=True)

            # Adding text to plot
            for i in range(0, len(lines)):
                vel_gal = c_in_km_per_s * ((1+lines[i].attrib["z_abs"])/(1+lines[i].attrib["z_gal"]) - 1)
                xticks_list = np.arange(-700, 800, 100)
                fig.axes[i].axvline(x=vel_gal, color='r', linestyle='--', linewidth=2)
                fig.axes[i].axhline(y=0, color='k', linestyle=(0, (1, 10)), linewidth=0.8)
                fig.axes[i].axhline(y=1, color='k', linestyle=(0, (1, 10)), linewidth=0.8)
                fig.axes[i].set_xticks(xticks_list)
                if len(lines) == 2:
                    fig.axes[i].text(-600, 0.24, "bval={:.4f}".format(lines[i].attrib["b"]))
                    fig.axes[i].text(-600, 0.17, "col={:.4f}".format(lines[i].attrib["logN"]))
                    fig.axes[i].text(-600, 0.1, "vel_gal={:.4f}".format(vel_gal))
                    fig.axes[i].text(600, 0.1, "z_abs={:.5f}".format(lines[i].attrib["z_abs"]))
                elif len(lines) == 3:
                    fig.axes[i].text(-600, 0.3, "bval={:.4f}".format(lines[i].attrib["b"]))
                    fig.axes[i].text(-600, 0.2, "col={:.4f}".format(lines[i].attrib["logN"]))
                    fig.axes[i].text(-600, 0.1, "vel_gal={:.4f}".format(vel_gal))
                    fig.axes[i].text(600, 0.1, "z_abs={:.5f}".format(lines[i].attrib["z_abs"]))
                elif len(lines) == 4:
                    fig.axes[i].text(-600, 0.38, "bval={:.4f}".format(lines[i].attrib["b"]))
                    fig.axes[i].text(-600, 0.26, "col={:.4f}".format(lines[i].attrib["logN"]))
                    fig.axes[i].text(-600, 0.14, "vel_gal={:.4f}".format(vel_gal))
                    fig.axes[i].text(600, 0.14, "z_abs={:.5f}".format(lines[i].attrib["z_abs"]))
                elif len(lines) == 5:
                    fig.axes[i].text(-600, 0.44, "bval={:.4f}".format(lines[i].attrib["b"]))
                    fig.axes[i].text(-600, 0.32, "col={:.4f}".format(lines[i].attrib["logN"]))
                    fig.axes[i].text(-600, 0.2, "vel_gal={:.4f}".format(vel_gal))
                    fig.axes[i].text(600, 0.2, "z_abs={:.5f}".format(lines[i].attrib["z_abs"]))
                else:
                    fig.axes[i].text(-600, 0.5, "bval={:.4f}".format(lines[i].attrib["b"]))
                    fig.axes[i].text(-600, 0.35, "col={:.4f}".format(lines[i].attrib["logN"]))
                    fig.axes[i].text(-600, 0.2, "vel_gal={:.4f}".format(vel_gal))
                    fig.axes[i].text(450, 0.1, "z_abs={:.5f}".format(lines[i].attrib["z_abs"]))
            ion_string = ''
            if aligned_df.trans.str.contains(low1).any() and aligned_df.trans.str.contains(low2).any():
                ion_string = f"{low1} and {low2}"
            elif aligned_df.trans.str.contains(low1).any():
                ion_string = low1
            elif aligned_df.trans.str.contains(low2).any():
                ion_string = low2

            if ion_string == '':
                string = f'{high} has no low-ion within range'
            else:
                string = f'{high} has low-ion(s) within range: '
            galaxy_row = galaxy_df.loc[galaxy_df["OBJECT"]==qso].loc[round(galaxy_df["z"], 7)==round(z_galaxy, 7)] # Make sure galaxy's redshift has same number of decimals included as z_galaxy
            rho_impact = galaxy_row.iloc[0]["rho_impact"]
            rho_rvir = galaxy_row.iloc[0]["rho_rvir"]
            mstars = galaxy_row.iloc[0]["mstars"]
            gal_type = galaxy_row.iloc[0]["gal_type"]

            title_type = ''
            if '/no-lows' in subdirs:
                title_type = 'No-Low'
            elif '/broads' in subdirs:
                title_type = 'Broad'
            elif '/narrows' in subdirs:
                title_type = 'Narrow'
            if len(lines) in (2, 3, 4, 5, 6):
                titlesize = 20
            else:
                titlesize = 14
            fig.axes[0].set_title(f'{string}{ion_string}  |  {title_type}\nQSO: {sightline}  |  z_gal={z_galaxy:.5f}  |  rho_impact={rho_impact:.4f}  |  rho_rvir={rho_rvir:.4f}  |  mstars={mstars:.2f}  |  gal_type={gal_type}', fontweight = 'bold', fontsize = titlesize, pad = 20)

            # Saving to pdf
            filename = f'{user_dir}stack_plots{subdirs[index]}/{sightline}_{z_galaxy:.5f}_stack_plot.pdf'
            fig.savefig(filename)
            plt.close(fig)


            # Saving stack plot data to csv
            z_gal_list = []
            z_abs_list = []
            vel_gal_list = []
            bval_list = []
            col_list = []
            transition_list = []
            for i in range(0, len(lines)):
                z_gal_list.append(lines[i].attrib["z_gal"])
                z_abs_list.append(lines[i].attrib["z_abs"])
                vel_gal = c_in_km_per_s * ((1+lines[i].attrib["z_abs"])/(1+lines[i].attrib["z_gal"]) - 1)
                vel_gal_list.append(vel_gal)
                bval_list.append(lines[i].attrib["b"])
                col_list.append(lines[i].attrib["logN"])
                transition_list.append(lines[i].attrib["trans"])
            alignment_type = ''
            if '/no-lows' in subdirs:
                alignment_type = 'No-Low'
            elif '/broads' in subdirs:
                alignment_type = 'Broad'
            elif '/narrows' in subdirs:
                alignment_type = 'Narrow'
            else:
                alignment_type = ''
            dat = {"filename" : filename, "QSO" : qso, "z_gal" : z_gal_list, "rho_impact" : rho_impact, "rho_rvir" : rho_rvir, "mstars" : mstars, "gal_type" : gal_type, "z_abs" : z_abs_list, "vel_gal" : vel_gal_list, "bval" : bval_list, "col" : col_list, "alignment_type" : alignment_type, "trans" : transition_list}
            stackplot_df = pd.DataFrame(data = dat)
            stackplot_df.to_csv(f'{user_dir}stack_plots_data/all_stack_plots{subdirs[index]}/{sightline}_{z_galaxy:.5f}_stack_plot_data.csv', sep=',')
    
    
    # Creating directories:
    directory = high
    path = os.path.join(user_dir, directory)
    if not os.path.exists(path):
        os.mkdir(path)
        print(f'\nCreated directory {directory}.\n')
    else:
        print(f'\n{directory} directory already exists. Overwriting existing {directory} directory.\n')
    
    user_dir = user_dir + high + '/' # Change in user_dir!
    
    directory = 'stack_plots'
    path = os.path.join(user_dir, directory)
    if not os.path.exists(path):
        os.mkdir(path)
    directory = 'stack_plots_data'
    path = os.path.join(user_dir, directory)
    if not os.path.exists(path):
        os.mkdir(path)
    directory = 'stack_plots_data/all_stack_plots'
    path = os.path.join(user_dir, directory)
    if not os.path.exists(path):
        os.mkdir(path)
    directory = 'stats'
    path = os.path.join(user_dir, directory)
    if not os.path.exists(path):
        os.mkdir(path)
    
    directory = 'stack_plots/all_alignments'
    path = os.path.join(user_dir, directory)
    if not os.path.exists(path):
        os.mkdir(path)
    directory = 'stack_plots/no-lows'
    path = os.path.join(user_dir, directory)
    if not os.path.exists(path):
        os.mkdir(path)
    directory = 'stack_plots/narrows'
    path = os.path.join(user_dir, directory)
    if not os.path.exists(path):
        os.mkdir(path)
    directory = 'stack_plots/broads'
    path = os.path.join(user_dir, directory)
    if not os.path.exists(path):
        os.mkdir(path)
#     directory = 'stack_plots/no-lows_strict'
#     path = os.path.join(user_dir, directory)
#     if not os.path.exists(path):
#         os.mkdir(path)
    
    directory = 'stack_plots_data/all_stack_plots/all_alignments'
    path = os.path.join(user_dir, directory)
    if not os.path.exists(path):
        os.mkdir(path)
    directory = 'stack_plots_data/all_stack_plots/no-lows'
    path = os.path.join(user_dir, directory)
    if not os.path.exists(path):
        os.mkdir(path)
    directory = 'stack_plots_data/all_stack_plots/narrows'
    path = os.path.join(user_dir, directory)
    if not os.path.exists(path):
        os.mkdir(path)
    directory = 'stack_plots_data/all_stack_plots/broads'
    path = os.path.join(user_dir, directory)
    if not os.path.exists(path):
        os.mkdir(path)
    
    
    # Preparing dataframe containing all high and low-ion transitions from finished all_csv that are within 500km/s of any galaxies in the galaxy_csv
    print(f'Identifying absorbers within 500km/s of galaxies...')
    galaxy_df = pd.read_csv(f'{galaxy_csv}', sep=',')
    all_df = pd.read_csv(f'{all_csv}', sep=',')
    all_df = all_df.sort_values(['QSO', 'z_comp'], ascending=True).reset_index(drop=True)

    vel_diff = 500
    # Criteria specific to H I high-ion:
    if high == 'H I':
        highall_df = all_df[all_df["trans"] == high]
        high_df = pd.DataFrame()
        for index, row in galaxy_df.iterrows():
            rho_impact_gal = row["rho_impact"]
            if rho_impact_gal <= 500: # rho_impact <= 500 kpc
                z_gal = row["z"]
                qso_gal = row["OBJECT"]
                for index, row in highall_df.iterrows():
                    col_abs = row["col"]
                    if col_abs > 14: # Hcol > 10^14
                        z_abs = row["zsys"]
                        qso_abs = row["QSO"]
                        z_diff = vel_diff*(1+z_gal)/c_in_km_per_s
                        if ((z_gal-z_diff) < z_abs < (z_gal+z_diff)) and (qso_abs == qso_gal):
                            row["z_gal"] = z_gal
                            high_df = high_df.append(row, ignore_index=True)
    # All other high-ions:
    else:
        highall_df = all_df[all_df["trans"] == high]
        high_df = pd.DataFrame()
        for index, row in galaxy_df.iterrows():
            z_gal = row["z"]
            qso_gal = row["OBJECT"]
            for index, row in highall_df.iterrows():
                z_abs = row["zsys"]
                qso_abs = row["QSO"]
                z_diff = vel_diff*(1+z_gal)/c_in_km_per_s
                if ((z_gal-z_diff) < z_abs < (z_gal+z_diff)) and (qso_abs == qso_gal):
                    row["z_gal"] = z_gal
                    high_df = high_df.append(row, ignore_index=True)

    low1all_df = all_df[all_df["trans"] == low1]
    low1_df = pd.DataFrame()
    for index, row in galaxy_df.iterrows():
        z_gal = row["z"]
        qso_gal = row["OBJECT"]
        for index, row in low1all_df.iterrows():
            z_abs = row["zsys"]
            qso_abs = row["QSO"]
            z_diff = vel_diff*(1+z_gal)/c_in_km_per_s
            if ((z_gal-z_diff) < z_abs < (z_gal+z_diff)) and (qso_abs == qso_gal):
                row["z_gal"] = z_gal
                low1_df = low1_df.append(row, ignore_index=True)
                
    low2all_df = all_df[all_df["trans"] == low2]
    low2_df = pd.DataFrame()
    for index, row in galaxy_df.iterrows():
        z_gal = row["z"]
        qso_gal = row["OBJECT"]
        for index, row in low2all_df.iterrows():
            z_abs = row["zsys"]
            qso_abs = row["QSO"]
            z_diff = vel_diff*(1+z_gal)/c_in_km_per_s
            if ((z_gal-z_diff) < z_abs < (z_gal+z_diff)) and (qso_abs == qso_gal):
                row["z_gal"] = z_gal
                low2_df = low2_df.append(row, ignore_index=True)
                
    dfs = [high_df, low1_df, low2_df]
    abs_df = pd.concat(dfs)
    abs_df = (
    abs_df
    .set_index(["QSO", "z_gal"])
    .sort_values(["QSO", "z_gal"], ascending=True)
    )
    print('Completed.\n')
          
    # Finding alignments:
    # Creating dataframe with all relevant ions in the qso
    print('Finding alignments and creating stack plots...')
    for qso, qso_df in abs_df.groupby(level=0):
        print(qso)
        if qso_df.trans.str.contains(high).any():
            high_df = (
                qso_df
                .query(f'trans == "%s"' % high)
                .sort_values("z_gal")
            )
            low1_df = (
                qso_df
                .query(f'trans == "%s"' % low1)
                .sort_values("z_gal")
            )
            low2_df = (
                qso_df
                .query(f'trans == "%s"' % low2)
                .sort_values("z_gal")
            )


            # Identifying all alignments between OVI and the low ions
            for z_high, z_highdf in high_df.groupby(level=1):
                #if math.floor(z_Odf.restwave[0]) == 1031: # Only plotting if OVI is 1031 and not 1037. This way we avoid printing two of the same plot
                #if z_highdf.index[0][1] == z_highdf.index[1][1]:
                aligned_df = pd.DataFrame()
                z_highdf["z_gal"] = z_high
                aligned_df = aligned_df.append(z_highdf)
                print(f'{high} {z_high}')

                if low1_df.index.size != 0:
                    for z_low1, z_low1df in low1_df.groupby(level=1):
                        if z_low1 == z_high:
                            z_low1df["z_gal"] = z_low1
                            aligned_df = aligned_df.append(z_low1df)
                            print(f'{low1} {z_low1}')

                if low2_df.index.size != 0:
                    for z_low2, z_low2df in low2_df.groupby(level=1):
                        if z_low2 == z_high:
                            z_low2df["z_gal"] = z_low2
                            aligned_df = aligned_df.append(z_low2df)
                            print(f'{low2} {z_low2}')


# No-Low:
                if (aligned_df.index.size != 0) and (aligned_df.trans.str.contains(low1).any() or aligned_df.trans.str.contains(low2).any()):
                    lines = []
                    indices = np.arange(0, aligned_df.index.size)
                    for ind in indices:
                        absline_df = aligned_df.iloc[ind]
                        absline_df["QSO"] = qso
                        line = make_absline(absline_df)
                        lines.append(line)
                    highvel_list = []
                    othervel_list = []
                    for line in lines:
                        if line.attrib['trans'] == high:
                            vel_gal = c_in_km_per_s * ((1+line.attrib["z_abs"])/(1+line.attrib["z_gal"]) - 1)
                            highvel_list.append(vel_gal)
                        else:
                            vel_gal = c_in_km_per_s * ((1+line.attrib["z_abs"])/(1+line.attrib["z_gal"]) - 1)
                            othervel_list.append(vel_gal)
                    vel_diffs = []
                    for i in range(0, len(highvel_list)):
                        for j in range(0, len(othervel_list)):
                            diff = highvel_list[i] - othervel_list[j]
                            diff = abs(diff)
                            vel_diffs.append(diff)
                    min_diff = min(vel_diffs)

                    if (min_diff > diff_tol) and not ((high == "O VI") and (qso == 'J1241+5721') and ((round(z_high, 3) == 0.147) or (round(z_high, 5) == 0.14650))):
                        keyword = "no-low"
                        splotter_plotter(high, low1, low2, lines, aligned_df, galaxy_df, qso, z_high, keyword, user_dir, spectra_dir)
# Broad/Narrow:
                    else:
                        b_vals = [] # Contains bvals for OVI
                        for line in lines:
                            if line.attrib['trans'] == high:
                                b_val = line.attrib["b"]
                                b_vals.append(b_val)
                        max_bval = max(b_vals)
# Broad:
                        if max_bval > 40*u.km/u.s:
                            keyword = "broad"
                            splotter_plotter(high, low1, low2, lines, aligned_df, galaxy_df, qso, z_high, keyword, user_dir, spectra_dir)
# Narrow:
                        else:
                            keyword = "narrow"
                            splotter_plotter(high, low1, low2, lines, aligned_df, galaxy_df, qso, z_high, keyword, user_dir, spectra_dir)

# No-Low (2nd sequence):
                elif (aligned_df.index.size != 0) and not (aligned_df.trans.str.contains(low1).any() or aligned_df.trans.str.contains(low2).any()):
                    lines = []
                    indices = np.arange(0, aligned_df.index.size)
                    for ind in indices:
                        absline_df = aligned_df.iloc[ind]
                        absline_df["QSO"] = qso
                        line = make_absline(absline_df)
                        lines.append(line)

                    keyword = "no-low"
                    splotter_plotter(high, low1, low2, lines, aligned_df, galaxy_df, qso, z_high, keyword, user_dir, spectra_dir)

    # Saving stats:
    all_dir = f'{user_dir}stack_plots/all_alignments'
    nolow_dir = f'{user_dir}stack_plots/no-lows'
    broad_dir = f'{user_dir}stack_plots/broads'
    narrow_dir = f'{user_dir}stack_plots/narrows'

    all_list = glob.glob(all_dir + '/*')
    nolow_list = glob.glob(nolow_dir + '/*')
    broad_list = glob.glob(broad_dir + '/*')
    narrow_list = glob.glob(narrow_dir + '/*')

    all_num = len(all_list)
    nolow_num = len(nolow_list)
    broad_num = len(broad_list)
    narrow_num = len(narrow_list)

    nolow_frac = nolow_num/all_num
    broad_frac = broad_num/all_num
    narrow_frac = narrow_num/all_num

    # Plotting stats        
    fig, ax = plt.subplots(1, 1)
    plt.style.use('fast')
    fig.set_size_inches(8,6)
    fig.tight_layout(w_pad = 5.0)
    ax.bar([1, 2, 3], [broad_num, narrow_num, nolow_num], tick_label = ['Broads', 'Narrows', 'No-Lows'],
            width = 0.8, color = ['tab:blue', 'tab:green', 'tab:red'])
    ax.set_xlabel('Alignment type', fontweight='bold', fontsize='15')
    ax.set_ylabel('Number of alignments', fontweight='bold', fontsize='15')
    ax.set_title(f'Kinematic Alignment Stats | High: {high}  Low: {low1}, {low2}', fontweight='bold', fontsize='20')
    ax.grid(linestyle='--')
    ax.xaxis.grid(False)
    ax.yaxis.grid(True, which='minor', linestyle='--', linewidth = 0.2, color='gray')
    ax.yaxis.grid(True, which='major', linestyle='-', linewidth = 0.5, color='k')
    if max([nolow_num, broad_num, narrow_num]) < 80:
        ax.yaxis.set_minor_locator(MultipleLocator(1))
        ax.yaxis.set_major_locator(MultipleLocator(5))
    ax.yaxis.set_major_formatter('{x:.0f}')
    ax.text(0.05, 0.95, f'Broad: {broad_num} ({broad_frac*100:.2f}%)\nNarrow: {narrow_num} ({narrow_frac*100:.2f}%)\nNo-Low: {nolow_num} ({nolow_frac*100:.2f}%)\nTotal: {all_num}', transform=ax.transAxes, fontsize=14, verticalalignment='top')

    fig.savefig(f'{user_dir}stats/alignment_stats.pdf', bbox_inches='tight')
    
    # Printing message:
    print(f"\n\nSuccessfully created stack plots for all matching high- and low-ion pairs in the directory \n"
          f"{user_dir}stack_plots/\n"
          f"Stack plot data can be found in the directory {user_dir}stack_plots_data/ \n\n"
          f"Splotter found {all_num} alignments in total whereas there are {broad_num} broads, {narrow_num} narrows, and {nolow_num} no-lows.\n\n"
          f"Done.")<jupyter_output><empty_output><jupyter_text>Example Code<jupyter_code>splotter("O VI", "SiIII", "C III", '/path/to/user_dir/', '/path/to/spectra/', '/path/to/all_vp.csv', '/path/to/galaxy_data_es.csv')
# INSTRUCTIONS:

# Leave "O VI", "SiIII" and "C III" as they are.

# Change the /path/to/ variables to match your directory:
# On my computer,
# '/path/to/user_dir/' looks like:
# '/Users/eriksolhaug/Dropbox/CGM2/Analysis/ErikSAnalysis/'

# '/path/to/spectra/' looks like:
# '/Users/eriksolhaug/Dropbox/CGM2/Analysis/ErikSAnalysis/spectra/'

# '/path/to/all_vp.csv' looks like:
# '/Users/eriksolhaug/Dropbox/CGM2/Analysis/ErikSAnalysis/finished_csv/all_vp.csv'
# NOTE: This one is not a directory but a file name!

# '/path/to/galaxy_data_es.csv' looks like:
# '/Users/eriksolhaug/Dropbox/CGM2/Analysis/ErikSAnalysis/data/galaxy_data_es.csv'
# NOTE: This one is not a directory but a file name!

# You can definitely set all the directory paths equal to the user directory path (likely where you placed this notebook) if this is easier!

# If the code runs successfully, then you should see a new subdirectory 'O VI' appear in your user directory ('user_dir')
# and the 'O VI' subdirectory should have a bunch of new pdf and csv files in them.
# The pdfs are 'stack plots' showing alignment between absorbers in the QSOs we are studying!! :)<jupyter_output>read_sets: Using set file -- 
  /Users/eriksolhaug/linetools/linetools/lists/sets/llist_v1.3.ascii
Loading abundances from Asplund2009
Abundances are relative by number on a logarithmic scale with H=12
{'fraction_non_alphanumeric': 0.12890459907019633, 'fraction_numerical': 0.016252353325392863, 'mean_word_length': 2.1366594360086766, 'pattern_counts': {'":': 3, '<': 19, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '599996', 'license': 'MIT', 'n_tokens_mistral': 9207, 'n_tokens_neox': 8614, 'n_words': 1864, 'path': 'Splotter_Code.ipynb', 'repo_name': 'SOLHAUG13/cgm_emission', 'size': '26027'}