#!/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/>.
"""
Very simple script to analyze GENIE MC to estimate the expected event rate from
atmospheric neutrino events before and after different cuts. To run it simply
pass the filename of a GENIE ntuple file:

    $ ./analyze-genie-mc mc_atm_nu_genie_010000_0.gst.root

Note: The input files should be in the GENIE "ntuple" format which you can
create by using the gntpc utility. For example:

    $ gntpc -i /sno/output/d2o-leta/mc_atmospherics/atmospheric_neutrinos_genie/root/mc_atm_nu_genie_010000_0.root -f gst

will create a new file mc_atm_nu_genie_010000_0.gst.root in the current working
directory in the ntuple file format.
"""
from __future__ import print_function, division
import ROOT
import numpy as np

# on retina screens, the default plots are way too small
# by using Qt5 and setting QT_AUTO_SCREEN_SCALE_FACTOR=1
# Qt5 will scale everything using the dpi in ~/.Xresources
import matplotlib
matplotlib.use("Qt5Agg")

def tick_formatter(x, pos):
    if x < 1:
        return '%.1f' % x
    else:
        return '%.0f' % x

# FIXME: What is this for the salt phase?
NEUTRON_DETECTION_EFFICIENCY = 0.5

# FIXME: What is the index for D2O?
INDEX_HEAVY_WATER = 1.3

# fractional energy resolution
# from Richie's thesis page 134
ENERGY_RESOLUTION = 0.05

def pdg_code_to_string(pdg):
    A = int(("%010i" % event.tgt)[-4:-1])
    Z = int(("%010i" % event.tgt)[-7:-4])
    if Z == 1:
        atom = 'H'
    elif Z == 8:
        atom = 'O'
    elif Z == 6:
        atom = 'C'
    else:
        raise NotImplementedError("unknown atom %i" % Z)
    return '%i%s' % (A,atom)

def get_reaction(event):
    reactants = []
    products = []
    if event.neu == 12:
        reactants.append('ve')
    elif event.neu == -12:
        reactants.append('vebar')
    elif event.neu == 14:
        reactants.append('vu')
    elif event.neu == -14:
        reactants.append('vubar')
    elif event.neu == 16:
        reactants.append('vt')
    elif event.neu == -16:
        reactants.append('vtbar')

    if event.hitnuc == 2212:
        reactants.append('p')
    elif event.hitnuc == 2112:
        reactants.append('n')
    elif event.hitnuc == 0:
        reactants.append(pdg_code_to_string(event.tgt))
    else:
        print("unknown nucleon %i" % event.hitnuc)

    if event.cc:
        if event.neu == 12:
            products.append('e-')
        elif event.neu == -12:
            products.append('e+')
        elif event.neu == 14:
            products.append('u-')
        elif event.neu == -14:
            products.append('u+')
        elif event.neu == 16:
            products.append('t-')
        elif event.neu == -16:
            products.append('t+')
    elif event.nc:
        if event.neu == 12:
            products.append('ve')
        elif event.neu == -12:
            products.append('vebar')
        elif event.neu == 14:
            products.append('vu')
        elif event.neu == -14:
            products.append('vubar')
        elif event.neu == 16:
            products.append('vt')
        elif event.neu == -16:
            products.append('vtbar')
    else:
        products.append("???")

    for pdg in event.pdgf:
        if pdg == 2112:
            products.append('n')
        elif abs(pdg) == 11:
            # e- or e+
            if pdg == 11:
                products.append('e-')
            else:
                products.append('e+')
        elif pdg == 22:
            # gamma
            products.append('gamma')
        elif pdg == 111:
            # pi0
            products.append('pi0')
        elif abs(pdg) == 211:
            # pi+/-
            if pdg == 211:
                products.append('pi+')
            else:
                products.append('pi-')
        elif abs(pdg) == 311:
            if pdg == 311:
                products.append('K0')
            else:
                products.append('K0bar')
        elif abs(pdg) == 321:
            # K+/-
            if pdg == 321:
                products.append('K+')
            else:
                products.append('K-')
        elif abs(pdg) == 3222:
            products.append('Sigma+')
        elif abs(pdg) == 3112:
            products.append('Sigma-')
        elif abs(pdg) == 3122:
            products.append('Delta')
        elif pdg == 2212:
            products.append('p')
        elif int(("%010i" % abs(pdg))[0]) == 1:
            products.append(pdg_code_to_string(pdg))
        else:
            print("unknown pdg code %i" % pdg)

    return ' + '.join(reactants) + ' -> ' + ' + '.join(products)

if __name__ == '__main__':
    import argparse
    import matplotlib.pyplot as plt
    from collections import Counter

    parser = argparse.ArgumentParser("script to analyze GENIE 'ntuple' ROOT files")
    parser.add_argument("filenames", nargs='+', help="GENIE ROOT files")
    args = parser.parse_args()

    bins = np.logspace(-1,2,100)

    El = []
    total_neutrons = []
    total_neutrons_detected = []
    E = []
    KE = []
    r = []
    total_nrings = []
    total_e_like_rings = []
    total_u_like_rings = []
    reactions = Counter()
    for filename in args.filenames:
        print("analyzing %s" % filename)
        f = ROOT.TFile(filename)
        T = f.Get("gst")

        for event in T:
            neutrons = 0
            nrings = 0
            e_like_rings = 0
            u_like_rings = 0
            ke = 0

            if event.cc:
                if abs(event.neu) == 12:
                    e_like_rings = 1
                else:
                    u_like_rings = 1
                nrings = 1
                ke += event.El
            elif event.nc:
                pass
            else:
                print("event is not cc or nc!")
                continue

            for i, pdg in enumerate(event.pdgf):
                if pdg == 2112:
                    neutrons += 1
                elif abs(pdg) == 11:
                    # e- or e+
                    if event.Ef[i] > 0.1:
                        # for now assume we only count rings from electrons
                        # with > 100 MeV
                        nrings += 1
                        e_like_rings += 1
                    ke += event.Ef[i]
                elif pdg == 22:
                    # gamma
                    if event.Ef[i] > 0.1:
                        # for now assume we only count rings from gammas with >
                        # 100 MeV
                        nrings += 1
                        e_like_rings += 1
                    ke += event.Ef[i]
                elif pdg == 111:
                    # pi0
                    nrings += 1
                    e_like_rings += 1
                    ke += event.Ef[i]
                elif abs(pdg) == 211:
                    # pi+/-
                    # momentum of ith particle in hadronic system
                    p = np.sqrt(event.pxf[i]**2 + event.pyf[i]**2 + event.pzf[i]**2)
                    # velocity of ith particle (in units of c)
                    # FIXME: is energy total energy or kinetic energy?
                    v = p/event.Ef[i]
                    if v > 1/INDEX_HEAVY_WATER:
                        # if the pion is above threshold, we assume that it
                        # produces 2 muon like rings
                        nrings += 2
                        u_like_rings += 2
                    else:
                        # if the pion is not above threshold, we assume that it
                        # produces 1 muon like ring
                        nrings += 1
                        u_like_rings += 1
                    # FIXME: should actually be a beta distribution
                    p = np.sqrt(event.pxf[i]**2 + event.pyf[i]**2 + event.pzf[i]**2)
                    m = np.sqrt(event.Ef[i]**2 - p**2)
                    ke += event.Ef[i] - m
                elif abs(pdg) in [2212,3222,311,321,3122,3112]:
                    # momentum of ith particle in hadronic system
                    p = np.sqrt(event.pxf[i]**2 + event.pyf[i]**2 + event.pzf[i]**2)
                    # velocity of ith particle (in units of c)
                    # FIXME: is energy total energy or kinetic energy?
                    v = p/event.Ef[i]
                    if v > 1/INDEX_HEAVY_WATER:
                        # above cerenkov threshold
                        nrings += 1
                        u_like_rings += 1
                        m = np.sqrt(event.Ef[i]**2 - p**2)
                        ke += event.Ef[i] - m
                elif int(("%010i" % abs(pdg))[0]) == 1:
                    # usually just excited 16O atom which won't produce a lot
                    # of light
                    pass
                else:
                    print("unknown pdg code %i" % pdg)
            total_neutrons.append(neutrons)
            total_neutrons_detected.append(np.random.binomial(neutrons,NEUTRON_DETECTION_EFFICIENCY))
            total_nrings.append(nrings)
            total_e_like_rings.append(e_like_rings)
            total_u_like_rings.append(u_like_rings)
            El.append(event.El)
            E.append(event.calresp0)
            KE.append(ke + np.random.randn()*ke*ENERGY_RESOLUTION)
            r.append(np.sqrt(event.vtxx**2 + event.vtxy**2 + event.vtxz**2))

            if total_neutrons_detected[-1] == 0 and nrings >= 2 and ((e_like_rings == 0) or (u_like_rings == 0)):
                reactions.update([get_reaction(event)])

    total = sum(reactions.values())
    for reaction, count in reactions.most_common(10):
        print("%.0f%% %s" % (count*100/total, reaction))

    El = np.array(El)
    total_neutrons = np.array(total_neutrons)
    total_neutrons_detected = np.array(total_neutrons_detected)
    E = np.array(E)
    KE = np.array(KE)
    r = np.array(r)
    total_nrings = np.array(total_nrings)
    total_e_like_rings = np.array(total_e_like_rings)
    total_u_like_rings = np.array(total_u_like_rings)

    cut1 = (total_neutrons_detected == 0)
    cut2 = (total_neutrons_detected == 0) & (total_nrings >= 2)
    cut3 = (total_neutrons_detected == 0) & (total_nrings >= 2) & ((total_e_like_rings == 0) | (total_u_like_rings == 0))

    El1 = El[cut1]
    El2 = El[cut2]
    El3 = El[cut3]
    E1 = E[cut1]
    E2 = E[cut2]
    E3 = E[cut3]
    KE1 = KE[cut1]
    KE2 = KE[cut2]
    KE3 = KE[cut3]

    plt.figure()
    bincenters = (bins[1:] + bins[:-1])/2
    x = np.repeat(bins,2)
    El_hist, _ = np.histogram(El, bins=bins)
    total_events = El_hist.sum()
    # FIXME: this is just a rough estimate of how many events we expect in 3
    # years based on Richie's thesis which says "Over the 306.4 live days of
    # the D2O phase we expect a total of 192.4 events within the acrylic vessel
    # and 504.5 events within the PSUP.
    El_hist = El_hist*230/total_events
    y = np.concatenate([[0],np.repeat(El_hist,2),[0]])
    El1_hist, _ = np.histogram(El1, bins=bins)
    El1_hist = El1_hist*230/total_events
    y1 = np.concatenate([[0],np.repeat(El1_hist,2),[0]])
    El2_hist, _ = np.histogram(El2, bins=bins)
    El2_hist = El2_hist*230/total_events
    y2 = np.concatenate([[0],np.repeat(El2_hist,2),[0]])
    El3_hist, _ = np.histogram(El3, bins=bins)
    El3_hist = El3_hist*230/total_events
    y3 = np.concatenate([[0],np.repeat(El3_hist,2),[0]])
    plt.plot(x, y, label="All events")
    plt.step(x, y1, where='mid', label="n")
    plt.step(x, y2, where='mid', label="n + nrings")
    plt.step(x, y3, where='mid', label="n + nrings + same flavor")
    plt.xlabel("Energy (GeV)")
    plt.ylabel("Events/year")
    plt.gca().set_xscale("log")
    plt.gca().xaxis.set_major_formatter(matplotlib.ticker.FuncFormatter(tick_formatter))
    plt.xlim((0.02,bins[-1]))
    plt.ylim((0,None))
    plt.legend()
    plt.title("Primary Lepton Energy")

    plt.figure()
    KE_hist, _ = np.histogram(KE, bins=bins)
    KE_signal, _ = np.histogram(np.random.randn(1000)*1.0*ENERGY_RESOLUTION + 1.0, bins=bins)
    total_events = KE_hist.sum()
    # FIXME: this is just a rough estimate of how many events we expect in 3
    # years based on Richie's thesis which says "Over the 306.4 live days of
    # the D2O phase we expect a total of 192.4 events within the acrylic vessel
    # and 504.5 events within the PSUP.
    KE_hist = KE_hist*230/total_events
    y = np.concatenate([[0],np.repeat(KE_hist,2),[0]])
    KE1_hist, _ = np.histogram(KE1, bins=bins)
    KE1_hist = KE1_hist*230/total_events
    y1 = np.concatenate([[0],np.repeat(KE1_hist,2),[0]])
    KE2_hist, _ = np.histogram(KE2, bins=bins)
    KE2_hist = KE2_hist*230/total_events
    y2 = np.concatenate([[0],np.repeat(KE2_hist,2),[0]])
    KE3_hist, _ = np.histogram(KE3, bins=bins)
    KE3_hist = KE3_hist*230/total_events
    y3 = np.concatenate([[0],np.repeat(KE3_hist,2),[0]])
    KE_signal = KE_signal*10/np.sum(KE_signal)
    y4 = np.concatenate([[0],np.repeat(KE_signal,2),[0]])
    plt.plot(x, y, label="All events")
    plt.plot(x, y1, label="n")
    plt.plot(x, y2, label="n + nrings")
    plt.plot(x, y3, label="n + nrings + same flavor")
    plt.plot(x, y4, label="1 GeV signal")
    plt.xlabel("Energy (GeV)")
    plt.ylabel(r"Expected Event Rate (year$^{-1}$)")
    plt.gca().set_xscale("log")
    plt.gca().xaxis.set_major_formatter(matplotlib.ticker.FuncFormatter(tick_formatter))
    plt.xlim((0.02,bins[-1]))
    plt.ylim((0,None))
    plt.legend()
    plt.title("Approximate Visible Energy")

    plt.figure()
    plt.hist(r, bins=np.linspace(0,8,100), histtype='step')
    plt.xlabel("R (m)")
    plt.title("Radius of Events")

    plt.figure()
    plt.hist(total_neutrons, bins=np.arange(11)-0.5, histtype='step')
    plt.xlabel("Number of Neutrons")
    plt.title("Number of Neutrons")

    plt.figure()
    plt.hist(total_nrings, bins=np.arange(11)-0.5, histtype='step')
    plt.xlabel("Number of Rings")
    plt.title("Number of Rings (approximate)")
    plt.show()
='n360' href='#n360'>360</a>
<a id='n361' href='#n361'>361</a>
<a id='n362' href='#n362'>362</a>
<a id='n363' href='#n363'>363</a>
<a id='n364' href='#n364'>364</a>
<a id='n365' href='#n365'>365</a>
<a id='n366' href='#n366'>366</a>
<a id='n367' href='#n367'>367</a>
<a id='n368' href='#n368'>368</a>
<a id='n369' href='#n369'>369</a>
<a id='n370' href='#n370'>370</a>
<a id='n371' href='#n371'>371</a>
<a id='n372' href='#n372'>372</a>
<a id='n373' href='#n373'>373</a>
<a id='n374' href='#n374'>374</a>
<a id='n375' href='#n375'>375</a>
<a id='n376' href='#n376'>376</a>
<a id='n377' href='#n377'>377</a>
<a id='n378' href='#n378'>378</a>
<a id='n379' href='#n379'>379</a>
<a id='n380' href='#n380'>380</a>
<a id='n381' href='#n381'>381</a>
<a id='n382' href='#n382'>382</a>
<a id='n383' href='#n383'>383</a>
<a id='n384' href='#n384'>384</a>
<a id='n385' href='#n385'>385</a>
<a id='n386' href='#n386'>386</a>
<a id='n387' href='#n387'>387</a>
<a id='n388' href='#n388'>388</a>
<a id='n389' href='#n389'>389</a>
<a id='n390' href='#n390'>390</a>
<a id='n391' href='#n391'>391</a>
<a id='n392' href='#n392'>392</a>
<a id='n393' href='#n393'>393</a>
<a id='n394' href='#n394'>394</a>
<a id='n395' href='#n395'>395</a>
<a id='n396' href='#n396'>396</a>
<a id='n397' href='#n397'>397</a>
<a id='n398' href='#n398'>398</a>
<a id='n399' href='#n399'>399</a>
<a id='n400' href='#n400'>400</a>
<a id='n401' href='#n401'>401</a>
<a id='n402' href='#n402'>402</a>
<a id='n403' href='#n403'>403</a>
<a id='n404' href='#n404'>404</a>
<a id='n405' href='#n405'>405</a>
<a id='n406' href='#n406'>406</a>
<a id='n407' href='#n407'>407</a>
<a id='n408' href='#n408'>408</a>
<a id='n409' href='#n409'>409</a>
<a id='n410' href='#n410'>410</a>
<a id='n411' href='#n411'>411</a>
<a id='n412' href='#n412'>412</a>
<a id='n413' href='#n413'>413</a>
<a id='n414' href='#n414'>414</a>
<a id='n415' href='#n415'>415</a>
<a id='n416' href='#n416'>416</a>
<a id='n417' href='#n417'>417</a>
<a id='n418' href='#n418'>418</a>
<a id='n419' href='#n419'>419</a>
<a id='n420' href='#n420'>420</a>
<a id='n421' href='#n421'>421</a>
<a id='n422' href='#n422'>422</a>
<a id='n423' href='#n423'>423</a>
<a id='n424' href='#n424'>424</a>
<a id='n425' href='#n425'>425</a>
<a id='n426' href='#n426'>426</a>
<a id='n427' href='#n427'>427</a>
<a id='n428' href='#n428'>428</a>
<a id='n429' href='#n429'>429</a>
<a id='n430' href='#n430'>430</a>
<a id='n431' href='#n431'>431</a>
<a id='n432' href='#n432'>432</a>
<a id='n433' href='#n433'>433</a>
<a id='n434' href='#n434'>434</a>
<a id='n435' href='#n435'>435</a>
<a id='n436' href='#n436'>436</a>
<a id='n437' href='#n437'>437</a>
<a id='n438' href='#n438'>438</a>
<a id='n439' href='#n439'>439</a>
<a id='n440' href='#n440'>440</a>
<a id='n441' href='#n441'>441</a>
<a id='n442' href='#n442'>442</a>
<a id='n443' href='#n443'>443</a>
<a id='n444' href='#n444'>444</a>
<a id='n445' href='#n445'>445</a>
<a id='n446' href='#n446'>446</a>
<a id='n447' href='#n447'>447</a>
<a id='n448' href='#n448'>448</a>
<a id='n449' href='#n449'>449</a>
<a id='n450' href='#n450'>450</a>
<a id='n451' href='#n451'>451</a>
<a id='n452' href='#n452'>452</a>
<a id='n453' href='#n453'>453</a>
<a id='n454' href='#n454'>454</a>
<a id='n455' href='#n455'>455</a>
<a id='n456' href='#n456'>456</a>
<a id='n457' href='#n457'>457</a>
<a id='n458' href='#n458'>458</a>
<a id='n459' href='#n459'>459</a>
<a id='n460' href='#n460'>460</a>
<a id='n461' href='#n461'>461</a>
<a id='n462' href='#n462'>462</a>
<a id='n463' href='#n463'>463</a>
<a id='n464' href='#n464'>464</a>
<a id='n465' href='#n465'>465</a>
<a id='n466' href='#n466'>466</a>
<a id='n467' href='#n467'>467</a>
<a id='n468' href='#n468'>468</a>
<a id='n469' href='#n469'>469</a>
<a id='n470' href='#n470'>470</a>
<a id='n471' href='#n471'>471</a>
<a id='n472' href='#n472'>472</a>
<a id='n473' href='#n473'>473</a>
<a id='n474' href='#n474'>474</a>
<a id='n475' href='#n475'>475</a>
<a id='n476' href='#n476'>476</a>
<a id='n477' href='#n477'>477</a>
<a id='n478' href='#n478'>478</a>
<a id='n479' href='#n479'>479</a>
<a id='n480' href='#n480'>480</a>
<a id='n481' href='#n481'>481</a>
<a id='n482' href='#n482'>482</a>
<a id='n483' href='#n483'>483</a>
<a id='n484' href='#n484'>484</a>
<a id='n485' href='#n485'>485</a>
<a id='n486' href='#n486'>486</a>
<a id='n487' href='#n487'>487</a>
<a id='n488' href='#n488'>488</a>
<a id='n489' href='#n489'>489</a>
<a id='n490' href='#n490'>490</a>
<a id='n491' href='#n491'>491</a>
<a id='n492' href='#n492'>492</a>
<a id='n493' href='#n493'>493</a>
<a id='n494' href='#n494'>494</a>
<a id='n495' href='#n495'>495</a>
<a id='n496' href='#n496'>496</a>
<a id='n497' href='#n497'>497</a>
<a id='n498' href='#n498'>498</a>
<a id='n499' href='#n499'>499</a>
<a id='n500' href='#n500'>500</a>
<a id='n501' href='#n501'>501</a>
<a id='n502' href='#n502'>502</a>
<a id='n503' href='#n503'>503</a>
<a id='n504' href='#n504'>504</a>
<a id='n505' href='#n505'>505</a>
<a id='n506' href='#n506'>506</a>
<a id='n507' href='#n507'>507</a>
<a id='n508' href='#n508'>508</a>
<a id='n509' href='#n509'>509</a>
<a id='n510' href='#n510'>510</a>
<a id='n511' href='#n511'>511</a>
<a id='n512' href='#n512'>512</a>
<a id='n513' href='#n513'>513</a>
<a id='n514' href='#n514'>514</a>
<a id='n515' href='#n515'>515</a>
<a id='n516' href='#n516'>516</a>
<a id='n517' href='#n517'>517</a>
<a id='n518' href='#n518'>518</a>
<a id='n519' href='#n519'>519</a>
<a id='n520' href='#n520'>520</a>
<a id='n521' href='#n521'>521</a>
<a id='n522' href='#n522'>522</a>
<a id='n523' href='#n523'>523</a>
<a id='n524' href='#n524'>524</a>
<a id='n525' href='#n525'>525</a>
<a id='n526' href='#n526'>526</a>
<a id='n527' href='#n527'>527</a>
<a id='n528' href='#n528'>528</a>
<a id='n529' href='#n529'>529</a>
<a id='n530' href='#n530'>530</a>
<a id='n531' href='#n531'>531</a>
<a id='n532' href='#n532'>532</a>
<a id='n533' href='#n533'>533</a>
<a id='n534' href='#n534'>534</a>
<a id='n535' href='#n535'>535</a>
<a id='n536' href='#n536'>536</a>
<a id='n537' href='#n537'>537</a>
<a id='n538' href='#n538'>538</a>
<a id='n539' href='#n539'>539</a>
<a id='n540' href='#n540'>540</a>
<a id='n541' href='#n541'>541</a>
<a id='n542' href='#n542'>542</a>
<a id='n543' href='#n543'>543</a>
<a id='n544' href='#n544'>544</a>
<a id='n545' href='#n545'>545</a>
<a id='n546' href='#n546'>546</a>
<a id='n547' href='#n547'>547</a>
<a id='n548' href='#n548'>548</a>
<a id='n549' href='#n549'>549</a>
<a id='n550' href='#n550'>550</a>
<a id='n551' href='#n551'>551</a>
<a id='n552' href='#n552'>552</a>
<a id='n553' href='#n553'>553</a>
<a id='n554' href='#n554'>554</a>
<a id='n555' href='#n555'>555</a>
<a id='n556' href='#n556'>556</a>
<a id='n557' href='#n557'>557</a>
<a id='n558' href='#n558'>558</a>
<a id='n559' href='#n559'>559</a>
<a id='n560' href='#n560'>560</a>
<a id='n561' href='#n561'>561</a>
<a id='n562' href='#n562'>562</a>
<a id='n563' href='#n563'>563</a>
<a id='n564' href='#n564'>564</a>
<a id='n565' href='#n565'>565</a>
<a id='n566' href='#n566'>566</a>
<a id='n567' href='#n567'>567</a>
<a id='n568' href='#n568'>568</a>
<a id='n569' href='#n569'>569</a>
<a id='n570' href='#n570'>570</a>
<a id='n571' href='#n571'>571</a>
<a id='n572' href='#n572'>572</a>
<a id='n573' href='#n573'>573</a>
<a id='n574' href='#n574'>574</a>
<a id='n575' href='#n575'>575</a>
<a id='n576' href='#n576'>576</a>
<a id='n577' href='#n577'>577</a>
<a id='n578' href='#n578'>578</a>
<a id='n579' href='#n579'>579</a>
<a id='n580' href='#n580'>580</a>
<a id='n581' href='#n581'>581</a>
<a id='n582' href='#n582'>582</a>
<a id='n583' href='#n583'>583</a>
<a id='n584' href='#n584'>584</a>
<a id='n585' href='#n585'>585</a>
<a id='n586' href='#n586'>586</a>
<a id='n587' href='#n587'>587</a>
<a id='n588' href='#n588'>588</a>
<a id='n589' href='#n589'>589</a>
<a id='n590' href='#n590'>590</a>
<a id='n591' href='#n591'>591</a>
<a id='n592' href='#n592'>592</a>
<a id='n593' href='#n593'>593</a>
<a id='n594' href='#n594'>594</a>
<a id='n595' href='#n595'>595</a>
<a id='n596' href='#n596'>596</a>
<a id='n597' href='#n597'>597</a>
<a id='n598' href='#n598'>598</a>
<a id='n599' href='#n599'>599</a>
<a id='n600' href='#n600'>600</a>
<a id='n601' href='#n601'>601</a>
<a id='n602' href='#n602'>602</a>
<a id='n603' href='#n603'>603</a>
<a id='n604' href='#n604'>604</a>
<a id='n605' href='#n605'>605</a>
<a id='n606' href='#n606'>606</a>
<a id='n607' href='#n607'>607</a>
<a id='n608' href='#n608'>608</a>
<a id='n609' href='#n609'>609</a>
<a id='n610' href='#n610'>610</a>
<a id='n611' href='#n611'>611</a>
<a id='n612' href='#n612'>612</a>
<a id='n613' href='#n613'>613</a>
<a id='n614' href='#n614'>614</a>
<a id='n615' href='#n615'>615</a>
<a id='n616' href='#n616'>616</a>
<a id='n617' href='#n617'>617</a>
<a id='n618' href='#n618'>618</a>
<a id='n619' href='#n619'>619</a>
<a id='n620' href='#n620'>620</a>
<a id='n621' href='#n621'>621</a>
<a id='n622' href='#n622'>622</a>
<a id='n623' href='#n623'>623</a>
<a id='n624' href='#n624'>624</a>
<a id='n625' href='#n625'>625</a>
<a id='n626' href='#n626'>626</a>
<a id='n627' href='#n627'>627</a>
<a id='n628' href='#n628'>628</a>
<a id='n629' href='#n629'>629</a>
<a id='n630' href='#n630'>630</a>
<a id='n631' href='#n631'>631</a>
<a id='n632' href='#n632'>632</a>
<a id='n633' href='#n633'>633</a>
<a id='n634' href='#n634'>634</a>
<a id='n635' href='#n635'>635</a>
<a id='n636' href='#n636'>636</a>
<a id='n637' href='#n637'>637</a>
<a id='n638' href='#n638'>638</a>
<a id='n639' href='#n639'>639</a>
<a id='n640' href='#n640'>640</a>
<a id='n641' href='#n641'>641</a>
<a id='n642' href='#n642'>642</a>
<a id='n643' href='#n643'>643</a>
<a id='n644' href='#n644'>644</a>
<a id='n645' href='#n645'>645</a>
<a id='n646' href='#n646'>646</a>
<a id='n647' href='#n647'>647</a>
<a id='n648' href='#n648'>648</a>
<a id='n649' href='#n649'>649</a>
<a id='n650' href='#n650'>650</a>
<a id='n651' href='#n651'>651</a>
<a id='n652' href='#n652'>652</a>
<a id='n653' href='#n653'>653</a>
<a id='n654' href='#n654'>654</a>
<a id='n655' href='#n655'>655</a>
<a id='n656' href='#n656'>656</a>
<a id='n657' href='#n657'>657</a>
<a id='n658' href='#n658'>658</a>
<a id='n659' href='#n659'>659</a>
<a id='n660' href='#n660'>660</a>
<a id='n661' href='#n661'>661</a>
<a id='n662' href='#n662'>662</a>
<a id='n663' href='#n663'>663</a>
<a id='n664' href='#n664'>664</a>
<a id='n665' href='#n665'>665</a>
<a id='n666' href='#n666'>666</a>
<a id='n667' href='#n667'>667</a>
<a id='n668' href='#n668'>668</a>
<a id='n669' href='#n669'>669</a>
<a id='n670' href='#n670'>670</a>
<a id='n671' href='#n671'>671</a>
<a id='n672' href='#n672'>672</a>
<a id='n673' href='#n673'>673</a>
<a id='n674' href='#n674'>674</a>
<a id='n675' href='#n675'>675</a>
<a id='n676' href='#n676'>676</a>
<a id='n677' href='#n677'>677</a>
<a id='n678' href='#n678'>678</a>
<a id='n679' href='#n679'>679</a>
<a id='n680' href='#n680'>680</a>
<a id='n681' href='#n681'>681</a>
<a id='n682' href='#n682'>682</a>
<a id='n683' href='#n683'>683</a>
<a id='n684' href='#n684'>684</a>
<a id='n685' href='#n685'>685</a>
<a id='n686' href='#n686'>686</a>
<a id='n687' href='#n687'>687</a>
<a id='n688' href='#n688'>688</a>
<a id='n689' href='#n689'>689</a>
<a id='n690' href='#n690'>690</a>
<a id='n691' href='#n691'>691</a>
<a id='n692' href='#n692'>692</a>
<a id='n693' href='#n693'>693</a>
<a id='n694' href='#n694'>694</a>
<a id='n695' href='#n695'>695</a>
<a id='n696' href='#n696'>696</a>
<a id='n697' href='#n697'>697</a>
<a id='n698' href='#n698'>698</a>
<a id='n699' href='#n699'>699</a>
<a id='n700' href='#n700'>700</a>
<a id='n701' href='#n701'>701</a>
<a id='n702' href='#n702'>702</a>
<a id='n703' href='#n703'>703</a>
<a id='n704' href='#n704'>704</a>
<a id='n705' href='#n705'>705</a>
<a id='n706' href='#n706'>706</a>
<a id='n707' href='#n707'>707</a>
<a id='n708' href='#n708'>708</a>
<a id='n709' href='#n709'>709</a>
</pre></td>
<td class='lines'><pre><code>