#include <thread>
#include <boost/filesystem.hpp>
#include <boost/range/adaptor/transformed.hpp>

#include "Database.h"
#include "tools.h"
#include "Time.h"
#include "Configuration.h"

#ifdef HAVE_HIGHLIGHT
#include "srchilite/sourcehighlight.h"
#include "srchilite/langmap.h"
#endif

#ifdef HAVE_ROOT
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#endif

using namespace std;
using boost::adaptors::transformed;

namespace fs = boost::filesystem;

// ------------------------------- Binning ----------------------------------

struct Binning : std::set<double>
{
    bool equidist { false };

    Binning() { }
    Binning(const Binning &m) : std::set<double>(m), equidist(m.equidist) { }
    Binning(size_t cnt, double mn, double mx) { set(cnt, mn, mx); }
    void set(size_t cnt, double mn, double mx)
    {
        if (cnt==0)
            return;

        if (empty())
            equidist = true;

        const double min = ::min(mn, mx);
        const double max = ::max(mn, mx);

        const double stp = (max-min)/cnt;

        for (int i=0; i<=cnt; i++)
            emplace(min+i*stp);
    }

    void add(double val)
    {
        emplace(val);
        equidist = false;
    }

    Binning &operator+=(const vector<double> &v)
    {
        if (!v.empty())
        {
            insert(v.cbegin(), v.cend());
            equidist = false;
        }
        return *this;
    }

    Binning operator+(const vector<double> &v)
    {
        return Binning(*this) += v;
    }

    void add(const double &val)
    {
        emplace(val);
        equidist = false;
    }

    string list() const
    {
        return boost::join(*this | transformed([](double d) { return std::to_string(d); }), ",");
    }

    string str() const
    {
        if (!equidist)
            return list();

        return to_string(size()-1)+":"+to_string(*begin())+","+to_string(*rbegin());
    }

    vector<double> vec() const { return vector<double>(begin(), end()); }

    //double operator[](size_t i) const { return vec().at(i); }
};

std::istream &operator>>(std::istream &in, Binning &m)
{
    const istreambuf_iterator<char> eos;
    string txt(istreambuf_iterator<char>(in), eos);

    const boost::regex expr(
                            "([0-9]+)[ /,:;]+"
                            "([-+]?[0-9]*[.]?[0-9]+([eE][-+]?[0-9]+)?)[ /,:;]+"
                            "([-+]?[0-9]*[.]?[0-9]+([eE][-+]?[0-9]+)?)"
                           );
    boost::smatch match;
    if (!boost::regex_match(txt, match, expr))
        throw runtime_error("Could not evaluate binning: "+txt);

    m = Binning(stoi(match[1].str()), stof(match[2].str()), stof(match[4].str()));

    return in;
}

std::ostream &operator<<(std::ostream &out, const Binning &m)
{
    return out << m.str();
}

// ---------------------------- Configuration -------------------------------

void SetupConfiguration(Configuration &conf)
{
    po::options_description control("Calcsource options");
    control.add_options()
        ("uri,u",          var<string>()
#if BOOST_VERSION >= 104200
         ->required()
#endif
         , "Database link as in\n\tuser:password@server[:port]/database[?compress=0|1].")
        ("out,o",            var<string>(conf.GetName()), "")
        ("dry-run",          po_switch(),   "Only write the queries to the query.log file. Internally overwrites uri with an empty string.")
        ;

    po::options_description binnings("Binnings");
    binnings.add_options()
        ("theta",            var<Binning>(Binning(90, 0, 90)),   "")
        ("theta-bin",        vars<double>(),                     "")
        ("esim",             var<Binning>(Binning(15, 2, 5)),    "")
        ("esim-bin",         vars<double>(),                     "")
        ("eest",             var<Binning>(Binning(15, 2, 5)),    "")
        ("eest-bin",         vars<double>(),                     "")
        ;

    po::options_description queries("Queries");
    queries.add_options()
        ("analysis",    var<string>("analysis.sql"),   "")
        ("data",        var<string>("data.sql"),       "")
        ("simulation",  var<string>("simulation.sql"), "")
        ;

    po::options_description preparation("Preparation");
    preparation.add_options()
        ("source-key", var<uint16_t>(5),          "")
        ("selector",   vars<string>(),            "")
        ("estimator",  var<string>()->required(), "")
        ("spectrum",   var<string>()->required(), "")
        ("env.*",      var<string>(),             "")
        ;

    po::options_description debug("Debug options");
    debug.add_options()
        ("print-connection", po_switch(),       "Print database connection information")
        ("print-queries",    po_switch(),       "")
        ("verbose,v",        var<uint16_t>(1),  "Verbosity (0: quiet, 1: default, 2: more, 3, ...)")
        ;

    //po::positional_options_description p;
    //p.add("file", 1); // The 1st positional options (n=1)

    conf.AddOptions(control);
    conf.AddOptions(binnings);
    conf.AddOptions(queries);
    conf.AddOptions(preparation);
    conf.AddOptions(debug);
    //conf.SetArgumentPositions(p);
}

void PrintUsage()
{
    //78 .............................................................................
    cout <<
        "spectrum - Calculate a spectrum with classical algorithms\n"
        "\n\n"
        "Usage: spectrum [-u URI] [options]\n"
        "\n"
        ;
    cout << endl;
}



// ----------------------------- Indentation --------------------------------

class sindent : public std::streambuf
{
    std::streambuf *sbuf   { nullptr };
    std::ostream   *owner  { nullptr };
    int             lastch { 0       }; // start-of-line
    std::string     str;

protected:
    virtual int overflow(int ch)
    {
        if (lastch=='\n' && ch != '\n' )
            sbuf->sputn(str.data(), str.size());

        return sbuf->sputc(lastch = ch);
    }

public:
    explicit sindent(std::streambuf *dest, uint16_t indent=0)
        : sbuf(dest), str(indent, ' ')
    {
    }

    explicit sindent(std::ostream& dest, uint16_t indent=0)
        : sbuf(dest.rdbuf()), owner(&dest), str(indent, ' ')
    {
        owner->rdbuf(this);
    }

    virtual ~sindent()
    {
        if (owner)
            owner->rdbuf(sbuf);
    }

    void set(uint16_t w) { str.assign(w, ' '); }
};

struct indent
{
    uint16_t w;
    indent(uint16_t _w=0) : w(_w) { }
};

std::ostream &operator<<(std::ostream &out, indent m)
{
    sindent *buf=dynamic_cast<sindent*>(out.rdbuf());
    if (buf)
        buf->set(m.w);
    return out;
}

struct separator
{
    string s;
    uint16_t n { 0 };
    separator(string _s="") : s(_s) { };
    separator(char c='-', uint16_t _n=78) : s(_n, c), n(_n) { };
};

std::ostream &operator<<(std::ostream &out, separator m)
{
    if (m.n)
        out << m.s;
    else
    {
        const uint8_t l = (78-m.s.size())/2-3;
        out << string(l<1?1:l, '-') << "={ " << m.s << " }=" << string(l<1?1:l, '-');
        if (m.s.size()%2)
            out << '-';
    }
    return out;
}

// ------------------------------- MySQL++ ----------------------------------

bool ShowWarnings(Database &connection)
{
    if (!connection.connected())
        return true;

    try
    {
        const auto resw =
            connection.query("SHOW WARNINGS").store();

        for (size_t i=0; i<resw.num_rows(); i++)
        {
            const mysqlpp::Row &roww = resw[i];

            cout << "\033[31m";
            cout << roww["Level"] << '[' << roww["Code"] << "]: ";
            cout << roww["Message"] << "\033[0m" << '\n';
        }
        cout << endl;
        return true;

    }
    catch (const exception &e)
    {
        cerr << "\nSHOW WARNINGS\n\n";
        cerr << "SQL query failed:\n" << e.what() << '\n' <<endl;
        return false;
    }
}

size_t Dump(ostream &out, Database &connection, const string &table)
{
    if (!connection.connected())
        return 0;

    out << connection.query("SHOW CREATE TABLE `"+table+"`").store()[0][1] << ";\n\n";

    const mysqlpp::StoreQueryResult res = connection.query("SELECT * FROM `"+table+"`").store();

    const string fields = boost::join(res.fields() | transformed([](const mysqlpp::Field &r) { return "`"+string(r.name())+"`"; }), ", ");

    out << "INSERT INTO `" << table << "` ( " << fields << " ) VALUES\n";
    out << boost::join(res | transformed([](const mysqlpp::Row &row) { return "  ( "+boost::join(row | transformed([](const mysqlpp::String &s) { return string(s);}), ", ")+" )";}), ",\n") << ";";
    out << "\n" << endl;

    return res.num_rows();
}

void PrintQuery(const string &query)
{
#ifdef HAVE_HIGHLIGHT
    stringstream qstr(query);
    srchilite::SourceHighlight sourceHighlight("esc256.outlang");
    sourceHighlight.setStyleFile("esc256.style");
    sourceHighlight.highlight(qstr, cout, "sql.lang");
    cout << endl;
#else
    cout << query << endl;
#endif
}

void CreateBinning(Database &connection, ostream &qlog, const Binning &bins, const string &name)
{
    mysqlpp::Query query0(&connection);
    query0 <<
        "CREATE TEMPORARY TABLE Binning" << name << "\n"
        "(\n"
        "   bin SMALLINT UNSIGNED NOT NULL,\n"
        "   lo FLOAT NOT NULL,\n"
        "   hi FLOAT NOT NULL,\n"
        "   PRIMARY KEY (bin) USING HASH\n"
        ")";

    qlog << query0 << ";\n" << endl;
    if (connection.connected())
        query0.execute();

    //connection.query("ALTER TABLE Binning"+name+" AUTO_INCREMENT=0").execute();

    mysqlpp::Query query1(&connection);
    query1 <<
        "INSERT INTO Binning" << name << " ( bin, lo, hi ) VALUES\n";

    // FIXME: Case of 1 and 2 bins

    // Bin  0: is the underflow bin...
    // Bin  N: is the overflow  bin...
    // Bin -1: if argument is NULL

    const auto vec = bins.vec();
    for (int i=1; i<vec.size()-2; i++)
        query1 << "  ( " << i << ", " << vec[i-1] << ", " << vec[i] << " ),\n";
    query1 << "  ( " << vec.size()-2 << ", " << vec[vec.size()-2] << ", " << vec[vec.size()-1] << " )\n";

    qlog << query1 << ";\n" << endl;

    if (connection.connected())
        cout << query1.execute().info() << endl;
    ShowWarnings(connection);
}

// ----------------------------- ROOT Histogram -----------------------------

struct Histogram
{
    // A workaround to be able to set a default also in C++11
    template<typename T, T def>
    struct Value
    {
        T t { def };
        Value() = default;
        Value(const T &_t) : t(_t) { }
        operator T() const { return t; }
    };

    string  name;
    string  title;
    string  dir;
    Binning binningx;
    Binning binningy;
    string  table;
    string  x;
    string  y;
    string  v;
    string  err;
    string  axisx;
    string  axisy;
    string  axisz;
    Value<bool,true> stats;
};

void WriteHistogram(Database &connection, const Histogram &hist)
{
#ifdef HAVE_ROOT
    if (!connection.connected())
        return;

    mysqlpp::Query query(&connection);
    query << "SELECT `%0:x` AS X,%1:y `%2:v` AS V%3:err FROM `%4:table`";
    query.parse();

    query.template_defaults["table"] = hist.table.c_str();

    query.template_defaults["x"] = hist.x.c_str();
    query.template_defaults["v"] = hist.v.c_str();
    if (!hist.y.empty())
        query.template_defaults["y"]   = (" `"+hist.y+"` AS Y,").c_str();
    if (!hist.err.empty())
        query.template_defaults["err"] = (", `"+hist.err+"` AS E").c_str();

    // PrintQuery(query.str());

    const mysqlpp::StoreQueryResult res = query.store();

    const auto vecx = hist.binningx.vec();
    const auto vecy = hist.binningy.vec();

    TH1 *h = 0;

    if (hist.y.empty())
    {
        h = hist.binningx.equidist ?
            new TH1D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.front(), vecx.back()) :
            new TH1D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.data());
    }
    else
    {
        if (hist.binningy.equidist)
        {
            h = hist.binningx.equidist ?
                new TH2D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.front(), vecx.back(), vecy.size()-1, vecy.front(), vecy.back()) :
                new TH2D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.data(), vecy.size()-1, vecy.front(), vecy.back());
        }
        else
        {
            h = hist.binningx.equidist ?
                new TH2D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.front(), vecx.back(), vecy.size()-1, vecy.front(), vecy.back()) :
                new TH2D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.data(), vecy.size()-1, vecy.front(), vecy.back());
        }
    }

    h->SetXTitle(hist.axisx.c_str());
    h->SetYTitle(hist.axisy.c_str());
    h->SetZTitle(hist.axisz.c_str());

    h->SetMarkerStyle(kFullDotMedium);

    h->SetStats(hist.stats);

    for (auto ir=res.cbegin(); ir!=res.cend(); ir++)
    {
        const auto &row = *ir;

        const uint32_t x = row["X"];
        const double   v = row["V"];
        if (x>h->GetNbinsX()+1)
            continue;

        try
        {
            const uint32_t y = row["Y"];
            if (y>h->GetNbinsY()+1)
                continue;

            h->SetBinContent(x, y, v);

        }
        catch (const mysqlpp::BadFieldName &e)
        {
            h->SetBinContent(x, v);
            try
            {
                h->SetBinError(x, double(row["E"]));
            }
            catch (const mysqlpp::BadFieldName &e)
            {
            }
        }
    }

    gFile->cd();
    TDirectory *dir = gFile->GetDirectory(hist.dir.c_str());
    if (dir)
        dir->cd();
    else
    {
        gFile->mkdir(hist.dir.c_str());
        gFile->cd(hist.dir.c_str());
    }
    h->Write();
    delete h;
#endif
}

// -------------------------------- Resources -------------------------------

#define RESOURCE(TYPE,RC) ([]() { \
    extern const char _binary_##RC##_start[]; \
    extern const char _binary_##RC##_end[];   \
    return TYPE(_binary_##RC##_start, _binary_##RC##_end); \
})()

string CreateResource(Configuration &conf, const string id, const string def, const string resource)
{
    string file = conf.Get<string>(id);

    if (file!=def)
    {
        file = conf.GetPrefixedString(id);
        cout << "Using " << file << "." << endl;
        return file;
    }

    file = conf.GetPrefixedString(id);

    cout << "Searching " << file << "... ";

    if (fs::exists(file))
    {
        cout << "found." << endl;
        return file;
    }

    cout << "creating!" << endl;

    ofstream(file) << resource;

    return file;
}

// ================================== MAIN ==================================


int main(int argc, const char* argv[])
{

    Time start;

    Configuration conf(argv[0]);
    conf.SetPrintUsage(PrintUsage);
    SetupConfiguration(conf);

    if (!conf.DoParse(argc, argv))
        return 127;

    // ----------------------------- Evaluate options --------------------------

    const string   uri     = conf.Get<bool>("dry-run") ? "" : conf.Get<string>("uri");
    const string   out     = conf.Get<string>("out");
    const uint16_t verbose = conf.Get<uint16_t>("verbose");

    const bool print_connection = conf.Get<bool>("print-connection");
    const bool print_queries    = conf.Get<bool>("print-queries");

    Binning binning_theta = conf.Get<Binning>("theta")+conf.Vec<double>("theta-bin");
    Binning binning_esim  = conf.Get<Binning>("esim");
    Binning binning_eest  = conf.Get<Binning>("eest");

    cout << "B\"theta\": " << binning_theta.str() << endl;
    cout << "B\"esim\":  " << binning_esim.str()  << endl;
    cout << "B\"eest\":  " << binning_eest.str()  << endl;

    const uint16_t source_key = conf.Get<uint16_t>("source-key");
    const string   where      = boost::join(conf.Vec<string>("selector"), " AND\n      ");
    const string   estimator  = conf.Get<string>("estimator");
    const string   spectrum   = conf.Get<string>("spectrum");
    const auto     env        = conf.GetOptions<string>("env.");

    cout << "\n";
    const string analysis_sql   = CreateResource(conf, "analysis",   "analysis.sql",   RESOURCE(std::string, spectrum_analysis_sql));
    const string data_sql       = CreateResource(conf, "data",       "data.sql",       RESOURCE(std::string, spectrum_data_sql));
    const string simulation_sql = CreateResource(conf, "simulation", "simulation.sql", RESOURCE(std::string, spectrum_simulation_sql));
    cout << endl;

    const string strzd = binning_theta.list();
    const string stre  = binning_eest.list();
    const string strth = binning_esim.list();

    // -------------------------------------------------------------------------
    // Checking for database connection

    Database connection(uri); // Keep alive while fetching rows

    if (!uri.empty())
    {
        if (verbose>0)
        {
            cout << "Connecting to database...\n";
            cout << "Client Version: " << mysqlpp::Connection().client_version() << endl;
        }

        try
        {
            connection.connected();
        }
        catch (const exception &e)
        {
            cerr << "SQL connection failed: " << e.what() << endl;
            return 1;
        }

        if (verbose>0)
        {
            cout << "Server Version: " << connection.server_version() << endl;
            cout << "Connected to " << connection.ipc_info() << endl;
        }

        if (print_connection)
        {
            try
            {
                const auto res1 = connection.query("SHOW STATUS LIKE 'Compression'").store();
                cout << "Compression of database connection is " << string(res1[0][1]) << endl;

                const auto res2 = connection.query("SHOW STATUS LIKE 'Ssl_cipher'").store();
                cout << "Connection to databases is " << (string(res2[0][1]).empty()?"UNENCRYPTED":"ENCRYPTED ("+string(res2[0][1])+")") << endl;
            }
            catch (const exception &e)
            {
                cerr << "\nSHOW STATUS LIKE 'Compression'\n\n";
                cerr << "SQL query failed:\n" << e.what() << endl;
                return 9;
            }
        }
    }

    // -------------------------------------------------------------------------
    // Create log streams

    ofstream qlog(out+".query.sql");
    ofstream flog(connection.connected() ? out+".dump.sql" : "");
    ofstream mlog(connection.connected() ? out+".C" : "");

    cout << "\n";
    cout << "Queries    will be logged  to " << out << ".query.sql\n";
    if (connection.connected())
    {
        cout << "Tables     will be dumped  to " << out << ".dump.sql\n";
        cout << "ROOT macro will be written to " << out << ".C\n";
    }

#ifdef HAVE_ROOT
    TFile root(connection.connected() ? (out+".hist.root").c_str() : "", "RECREATE");
    if (connection.connected())
    {
        if (root.IsZombie())
            return 10;
        cout << "Histograms will be written to " << out << ".hist.root\n";
    }
#endif

    cout << endl;

    // FIMXE: Implement SYNTAX check on spectrum, estimator and selector

    // -------------------------------------------------------------------
    // --------------------------- 0: Binnings ---------------------------

    cout << separator("Binnings") << '\n';

    CreateBinning(connection, qlog, binning_theta, "Theta");
    CreateBinning(connection, qlog, binning_eest,  "EnergyEst");
    CreateBinning(connection, qlog, binning_esim,  "EnergySim");

    Dump(flog, connection, "BinningTheta");
    Dump(flog, connection, "BinningEnergyEst");
    Dump(flog, connection, "BinningEnergySim");

    // -------------------------------------------------------------------
    // --------------------------- 1: DataFiles --------------------------
    cout << separator("DataFiles") << '\n';

    Time start1;

    /* 01:get-data-files.sql */
    mysqlpp::Query query1(&connection);
    query1 <<
        "CREATE TEMPORARY TABLE DataFiles\n"
        "(\n"
        "   FileId INT UNSIGNED NOT NULL,\n"
        "   PRIMARY KEY (FileId) USING HASH\n"
        ") ENGINE=Memory\n"
        "AS\n"
        "(\n"
        "   SELECT\n"
        "      FileId\n"
        "   FROM\n"
        "      factdata.RunInfo\n"
        "   WHERE\n"
        "      fRunTypeKEY=1 AND fSourceKEY=%100:source AND\n"
        "      %101:where\n"
        "   ORDER BY\n"
        "      FileId\n"  // In order: faster
        ")";

    query1.parse();
    for (auto it=env.cbegin(); it!=env.cend(); it++)
        query1.template_defaults[it->first.c_str()] = it->second.c_str();

    query1.template_defaults["source"] = to_string(source_key).c_str();
    query1.template_defaults["where"]  = where.c_str();

    if (print_queries)
        PrintQuery(query1.str());

    qlog << query1 << ";\n" << endl;
    if (connection.connected())
    {
        cout << query1.execute().info() << endl;
        ShowWarnings(connection);
        Dump(flog, connection, "DataFiles");

        const auto sec1 = Time().UnixTime()-start1.UnixTime();
        cout << "Execution time: " << sec1 << "s\n" << endl;
    }

    // FIXME: Setup Zd binning depending on Data

    // -------------------------------------------------------------------
    // ------------------------ 2: ObservationTime -----------------------

    cout << separator("ObservationTime") << '\n';

    Time start2;

    /* 02:get-observation-time.sql */
    mysqlpp::Query query2(&connection);
    query2 <<
        "CREATE TEMPORARY TABLE ObservationTime\n"
        "(\n"
        "   `.theta` SMALLINT UNSIGNED NOT NULL,\n"
        "   OnTime FLOAT NOT NULL,\n"
        "   PRIMARY KEY (`.theta`) USING HASH\n"
        ") ENGINE=Memory\n"
        "AS\n"
        "(\n"
        "   SELECT\n"
        "      INTERVAL(fZenithDistanceMean, %100:bins) AS `.theta`,\n"
        "      SUM(TIME_TO_SEC(TIMEDIFF(fRunStop,fRunStart))*fEffectiveOn) AS OnTime\n"
        "   FROM\n"
        "      DataFiles\n"
        "   LEFT JOIN \n"
        "      factdata.RunInfo USING (FileId)\n"
        "   GROUP BY\n"
        "      `.theta`\n"
        "   ORDER BY\n"
        "      `.theta`\n"
        ")";

    query2.parse();
    for (auto it=env.cbegin(); it!=env.cend(); it++)
        query2.template_defaults[it->first.c_str()] = it->second.c_str();

    query2.template_defaults["bins"] = strzd.c_str();

    if (print_queries)
        PrintQuery(query2.str());

    qlog << query2 << ";\n" << endl;
    if (connection.connected())
    {
        cout << query2.execute().info() << endl;
        ShowWarnings(connection);
        Dump(flog, connection, "ObservationTime");

        const auto sec2 = Time().UnixTime()-start2.UnixTime();
        cout << "Execution time: " << sec2 << "s\n\n";
    }

    /*

    SELECT
       min(first_number)   as first_number,
       max(last_number)    as last_number
    from
    (
       select
          first_number,
          last_number
          sum(grp) over(order by first_number) as grp
       from
       (
          select
              bin AS first_number,
              LEAD(bin) AS last_number
              IF(first_number <> lag(bin, 2) over (order by bin) + 1, 1, 0) as grp
          from
              t1
       )
    )
    group by grp
    order by 1
    */

    /*
WITH Hallo AS (WITH MyTest AS (WITH MyTable AS (SELECT bin FROM Test
WHERE cnt>0 GROUP BY bin ORDER BY bin) SELECT bin AS first_number,
LAG(bin) OVER(ORDER BY bin) AS prev, LEAD(bin) OVER(ORDER BY bin) AS
last_number, IF(bin <> (LAG(bin) OVER(ORDER BY bin)+1), 1, 0) AS
grp_prev, IF(bin <> (LEAD(bin) OVER(ORDER BY bin)-1), 1, 0) as grp_nxt
FROM MyTable ) SELECT first_number, last_number, sum(grp_nxt)
over(order by first_number) as grp FROM MyTest) SELECT * FROM Hallo


    WITH MyTable AS
(SELECT bin FROM Test WHERE cnt>0 GROUP BY bin ORDER BY bin)
SELECT bin, LAG(bin) OVER(ORDER BY bin) AS prev, LEAD(bin) OVER(ORDER BY bin) AS nxt, 
IF(bin = (LAG(bin) OVER(ORDER BY bin)+1), 1, 0) AS grp_prev, IF(bin = (LEAD(bin) OVER(ORDER BY bin)-1), 1, 0) as grp_nxt  FROM MyTable

"-1","0"
"0","0"
"2","3"       2-4
"3","5"
"4","6"
"5","0"
"6","0"
"7","7"      7-7
"8","0"
"9","9"      9-12
"10","10"
"11","11"
"12","12"
"13","0"
"14","0"
"15","0"
"16","16"   16-19
"17","17"
"18","18"
"19","19"
    */

    // -------------------------------------------------------------------
    // ------------------------ 3: MonteCarloFiles -----------------------

    cout << separator("MonteCarloFiles") << '\n';

    Time start3;

    /* 02:get-monte-carlo-files.sql */
    mysqlpp::Query query3(&connection);
    query3 <<
        "CREATE TEMPORARY TABLE MonteCarloFiles\n"
        "(\n"
        "   FileId INT UNSIGNED NOT NULL,\n"
        "   PRIMARY KEY (FileId) USING HASH\n"
        ") ENGINE=Memory\n"
        "AS\n"
        "(\n"
        "   SELECT\n"
        "      FileId\n"
        "   FROM\n"
        "      ObservationTime\n"
        "   LEFT JOIN\n"
        "      BinningTheta ON `.theta`=bin\n"
        "   LEFT JOIN\n"
        "      factmc.RunInfoMC\n"
        "   ON\n"
        //"      ThetaMin BETWEEN lo AND hi OR ThetaMax BETWEEN lo AND hi\n" // Includes BOTH edges
        "      (ThetaMin>=lo AND ThetaMin<hi) OR (ThetaMax>lo AND ThetaMax<=hi)\n"
        "   WHERE\n"
        "      PartId=1 AND\n"
        "      FileId%%2=0\n"
        "   ORDER BY\n"
        "      FileId\n" // In order: faster
        ")";
/*
        "   SELECT\n"
        "      FileId\n"
        "   FROM\n"
        "      factmc.RunInfoMC\n"
        "   WHERE\n"
        "      PartId=1 AND\n"
        "      FileId%%2=0 AND\n"
        "      (%0:range)\n"
        ")";
*/

    query3.parse();
    for (auto it=env.cbegin(); it!=env.cend(); it++)
        query3.template_defaults[it->first.c_str()] = it->second.c_str();

    if (print_queries)
        PrintQuery(query3.str());

    qlog << query3 << ";\n" << endl;
    if (connection.connected())
    {
        cout << query3.execute().info() << endl;
        ShowWarnings(connection);
        Dump(flog, connection, "MonteCarloFiles");

        const auto sec3 = Time().UnixTime()-start3.UnixTime();
        cout << "Execution time: " << sec3 << "s\n\n";
    }

    // -------------------------------------------------------------------
    // ----------------------- 3b: Monte Carlo Area ----------------------

    cout << separator("MonteCarloArea") << '\n';

    Time start3b;

    /* 02:get-monte-carlo-files.sql */
    mysqlpp::Query query3b(&connection);
    query3b <<
        "CREATE TEMPORARY TABLE MonteCarloArea ENGINE=Memory\n"
        "AS\n"
        "(\n"
        "   SELECT\n"
        "      MIN(`CSCAT[1]`) AS MinImpactLo,\n"
        "      MAX(`CSCAT[1]`) AS MaxImpactLo,\n"
        "      MIN(`CSCAT[2]`) AS MinImpactHi,\n"
        "      MAX(`CSCAT[2]`) AS MaxImpactHi\n"
        "   FROM\n"
        "      MonteCarloFiles\n"
        "   LEFT JOIN\n"
        "      factmc.CorsikaSetup ON FileId=RUNNR\n"
        "   GROUP BY\n"
        "      `CSCAT[1]`, `CSCAT[2]`\n"
        "   ORDER BY\n"
        "      MaxImpactHi\n"
        ")";

    query3b.parse();
    for (auto it=env.cbegin(); it!=env.cend(); it++)
        query3b.template_defaults[it->first.c_str()] = it->second.c_str();

    if (print_queries)
        PrintQuery(query3b.str());

    qlog << query3b << ";\n" << endl;
    if (connection.connected())
    {
        cout << query3b.execute().info() << endl;
        ShowWarnings(connection);
        if (Dump(flog, connection, "MonteCarloArea")!=1)
        {
            cerr << "Impact range inconsistent!" << endl;
            return 1;
        }

        const auto sec3b = Time().UnixTime()-start3b.UnixTime();
        cout << "Execution time: " << sec3b << "s\n\n";
    }


    // -------------------------------------------------------------------
    // -------------------------- 4: ThetaHist ---------------------------

    cout << separator("ThetaHist") << '\n';

    Time start4;

    /* 02:get-theta-distribution.sql */
    mysqlpp::Query query4(&connection);
    query4 <<
        "CREATE TEMPORARY TABLE ThetaHist\n"
        "(\n"
        "   `.theta`    SMALLINT UNSIGNED NOT NULL,\n"
        "   lo          DOUBLE            NOT NULL COMMENT 'Lower edge of zenith distance bin in degree',\n"
        "   hi          DOUBLE            NOT NULL COMMENT 'Upper edge of zenith distance bin in degree',\n"
        "   CountN      INT UNSIGNED      NOT NULL,\n"
        "   ErrCountN   DOUBLE            NOT NULL,\n"
        "   OnTime      FLOAT             NOT NULL,\n"
        "   ZdWeight    DOUBLE            NOT NULL COMMENT 'tau(delta theta)',\n"
        "   ErrZdWeight DOUBLE            NOT NULL COMMENT 'sigma(tau)',\n"
        "   PRIMARY KEY (`.theta`) USING HASH\n"
        ") ENGINE=Memory\n"
        "AS\n"
        "(\n"
        "   WITH EventCount AS\n"
        "   (\n"
        "      SELECT\n"
        "         INTERVAL(DEGREES(Theta), %100:bins) AS `.theta`,\n"
        "         COUNT(*) AS CountN\n"
        "      FROM\n"
        "         MonteCarloFiles\n"
        "      LEFT JOIN\n"
        "         factmc.OriginalMC USING(FileId)\n"
        "      GROUP BY\n"
        "         `.theta`\n"
        "   )\n"
        "   SELECT\n"
        "      `.theta`, lo, hi,\n"
        "      CountN,\n"
        "      SQRT(CountN) AS ErrCountN,\n"
        "      OnTime,\n"
        "      OnTime/CountN AS ZdWeight,\n"
        "      (OnTime/CountN)*SQRT(POW(1/300, 2) + 1/CountN) AS ErrZdWeight\n"
        "   FROM\n"
        "      ObservationTime\n"
        "   LEFT JOIN\n"
        "      EventCount USING(`.theta`)\n"
        "   LEFT JOIN\n"
        "      BinningTheta ON `.theta`=bin\n"
        "   ORDER BY\n"
        "      `.theta`\n"
        ")";

    query4.parse();
    for (auto it=env.cbegin(); it!=env.cend(); it++)
        query4.template_defaults[it->first.c_str()] = it->second.c_str();

    query4.template_defaults["bins"] = strzd.c_str();

    if (print_queries)
        PrintQuery(query4.str());

    qlog << query4 << ";\n" << endl;
    if (connection.connected())
    {
        cout << query4.execute().info() << endl;
        ShowWarnings(connection);
        Dump(flog, connection, "ThetaHist");

        const auto sec4 = Time().UnixTime()-start4.UnixTime();
        cout << "Execution time: " << sec4 << "s\n";


        const mysqlpp::StoreQueryResult res4 = connection.query("SELECT * FROM ThetaHist").store();

        cout << "  Bin   MC Counts     OnTime ZdWeight\n";
        const auto bins = binning_theta.vec();
        for (auto ir=res4.cbegin(); ir!=res4.cend(); ir++)
        {
            const mysqlpp::Row &row = *ir;

            const uint32_t bin = row[".theta"];
            cout << setw(5) << bins[bin] << ": " << setw(10) << row["CountN"] << " " << setw(10) << row["OnTime"] << " " << setw(10) << row["ZdWeight"] << '\n';
        }
        cout << endl;
    }

    WriteHistogram(connection, {
             .dir      = "Zd",
             .name     = "OnTime",
             .title    = "Effective on time",
             .binningx = binning_theta,
             .table    = "ThetaHist",
             .x        = ".theta",
             .v        = "OnTime",
             .axisx    = "Zenith Distance #theta [#circ]",
             .axisy    = "Eff. on time [s]"
         });

    WriteHistogram(connection, {
             .dir      = "Zd",
             .name     = "CountN",
             .title    = "Simulated Zenith Distance",
             .binningx = binning_theta,
             .table    = "ThetaHist",
             .x        = ".theta",
             .v        = "CountN",
             .err      = "ErrCountN",
             .axisx    = "Zenith Distance #theta [#circ]",
             .axisy    = "Counts"
         });

    WriteHistogram(connection, {
             .dir      = "Zd",
             .name     = "ZdWeight",
             .title    = "Zenith Distance Weight",
             .binningx = binning_theta,
             .table    = "ThetaHist",
             .x        = ".theta",
             .v        = "ZdWeight",
             .err      = "ErrZdWeight",
             .axisx    = "Zenith Distance #theta [#circ]",
             .axisy    = "Weight [s]"
         });


    // -------------------------------------------------------------------
    // ------------------------- 5: AnalysisData -------------------------

    cout << separator("AnalysisData") << '\n';

    Time start5;

    /* 02:analyze-data.sql */
    mysqlpp::Query query5(&connection);
    sindent indent5(query5);
    query5 <<
        "CREATE TEMPORARY TABLE AnalysisData\n"
        "(\n"
        "   `.energy`       SMALLINT UNSIGNED NOT NULL,\n"
        "   `Signal`        DOUBLE            NOT NULL,\n"
        "   `Background`    DOUBLE            NOT NULL,\n"
        "   `Excess`        DOUBLE            NOT NULL,\n"
        "   `Significance`  DOUBLE            NOT NULL,\n"
        "   `ErrExcess`     DOUBLE            NOT NULL,\n"
        "   PRIMARY KEY (`.energy`) USING HASH\n"
        ") ENGINE=Memory\n"
        "AS\n"
        "(\n"
        "   WITH Excess AS\n"
        "   (\n"                          << indent(6)
        << ifstream(analysis_sql).rdbuf() << indent(0) <<
        "   ),\n"                         <<
        "   Result AS\n"
        "   (\n"                          << indent(6)
        << ifstream(data_sql).rdbuf()     << indent(0) << // Must end with EOL and not in the middle of a comment
        "   )\n"
        "   SELECT * FROM Result\n"
        ")";

    query5.parse();
    for (auto it=env.cbegin(); it!=env.cend(); it++)
        query5.template_defaults[it->first.c_str()] = it->second.c_str();

    //query5.template_defaults["columns"]   = "FileId, EvtNumber,";
    query5.template_defaults["columns"]   = "";
    query5.template_defaults["files"]     = "DataFiles";
    query5.template_defaults["runinfo"]   = "factdata.RunInfo";
    query5.template_defaults["events"]    = "factdata.Images";
    query5.template_defaults["positions"] = "factdata.Position";

    query5.template_defaults["bins"]      = stre.c_str();
    query5.template_defaults["estimator"] = estimator.c_str();

    if (print_queries)
        PrintQuery(query5.str());

    qlog << query5 << ";\n" << endl;
    if (connection.connected())
    {
        cout << query5.execute().info() << endl;
        ShowWarnings(connection);
        Dump(flog, connection, "AnalysisData");

        const auto sec5 = Time().UnixTime()-start5.UnixTime();
        cout << "Execution time: " << sec5 << "s\n";


        const mysqlpp::StoreQueryResult res5 = connection.query("SELECT * FROM AnalysisData").store();

        cout << "       Bin     Signal   Background   Excess   Significance   Error" << endl;
        for (auto row=res5.cbegin(); row!=res5.cend(); row++)
        {
            for (auto it=row->begin(); it!=row->end(); it++)
                cout << setw(10) << *it << " ";
            cout << '\n';
        }
        cout << endl;
    }


    // -------------------------------------------------------------------
    // --------------------------- 6: ResultMC ---------------------------

    cout << separator("ResultMC") << '\n';

    Time start6;

    /* 02:analyze-simulation.sql */

    // This table combines the analysis results vs. Binning in Estimated Energy and Simulated Energy
    mysqlpp::Query query6(&connection);
    sindent indent6(query6);
    query6 <<
        "CREATE TEMPORARY TABLE AnalysisMC\n"
        "(\n"
        "   `.energyest`  SMALLINT UNSIGNED NOT NULL,\n"
        "   `.energysim`  SMALLINT UNSIGNED NOT NULL,\n"
        "   SignalW       DOUBLE            NOT NULL,\n"
        "   BackgroundW   DOUBLE            NOT NULL,\n"
        "   ThresholdW    DOUBLE            NOT NULL,\n"
        "   SignalN       DOUBLE            NOT NULL,\n"
        "   BackgroundN   DOUBLE            NOT NULL,\n"
        "   ThresholdN    DOUBLE            NOT NULL,\n"
        "   ExcessW       DOUBLE            NOT NULL,\n"
        "   ExcessN       DOUBLE            NOT NULL,\n"
        "   ErrExcessW    DOUBLE            NOT NULL,\n"
        "   ErrExcessN    DOUBLE            NOT NULL,\n"
        "   BiasEst       DOUBLE            NOT NULL,\n"
        "   BiasSim       DOUBLE            NOT NULL,\n"
        "   ResolutionEst DOUBLE,\n"
        "   ResolutionSim DOUBLE,\n"
        "   INDEX (`.energyest`) USING HASH,\n"
        "   INDEX (`.energysim`) USING HASH\n"
        ") ENGINE=Memory\n"
        "AS\n"
        "(\n"
        "   WITH Excess AS\n"
        "   (\n"                            << indent(6)
        << ifstream(analysis_sql).rdbuf()   << indent(0) <<
        "   ),\n"                           <<
        "   Result AS\n"
        "   (\n"                            << indent(6)
        << ifstream(simulation_sql).rdbuf() << indent(0) << // Must end with EOL and not in the middle of a comment
        "   )\n"
        "   SELECT * FROM Result\n"
        ")";

    query6.parse();
    for (auto it=env.cbegin(); it!=env.cend(); it++)
        query6.template_defaults[it->first.c_str()] = it->second.c_str();

    //query6.template_defaults["columns"]   = "FileId, EvtNumber, CorsikaNumReuse,";
    query6.template_defaults["columns"]   = "Zd, Energy, SpectralIndex,";
    query6.template_defaults["files"]     = "MonteCarloFiles";
    query6.template_defaults["runinfo"]   = "factmc.RunInfoMC";
    query6.template_defaults["events"]    = "factmc.EventsMC";
    query6.template_defaults["positions"] = "factmc.PositionMC";

    query6.template_defaults["energyest"] = stre.c_str();
    query6.template_defaults["energysim"] = strth.c_str();
    query6.template_defaults["theta"]     = strzd.c_str();
    query6.template_defaults["spectrum"]  = spectrum.c_str();
    query6.template_defaults["estimator"] = estimator.c_str();

    if (print_queries)
        PrintQuery(query6.str());

    qlog << query6 << ";\n" << endl;
    if (connection.connected())
    {
        cout << query6.execute().info() << endl;
        ShowWarnings(connection);
        Dump(flog, connection, "AnalysisMC");

        const auto sec6 = Time().UnixTime()-start6.UnixTime();
        cout << "Execution time: " << sec6 << "s\n\n";
    }


    // -------------------------------------------------------------------
    // ----------------------- 7: SimulatedSpectrum ----------------------

    cout << separator("SimulatedSpectrum") << '\n';

    // [lo^(1+gammaa) - hi^(1+gamma)] / (1+gamma)

    Time start7;
    /* 02:get-corsika-events.sql */

    // FIXME: Theta weights?
    // FIXME: energysim binning
    mysqlpp::Query query7(&connection);
    query7 <<
        "CREATE TEMPORARY TABLE SimulatedSpectrum\n"
        "(\n"
        "   `.energy` SMALLINT UNSIGNED NOT NULL COMMENT 'Bin Index [MC Energy]',\n"
        "   CountN    DOUBLE            NOT NULL,\n"  // COMMENT 'Number of events',\n"
        "   CountW    DOUBLE            NOT NULL,\n"  // COMMENT 'Sum of weights, reweighted sum',\n"
        "   CountW2   DOUBLE            NOT NULL,\n"  // COMMENT 'Sum of square of weights'\n"
        "   PRIMARY KEY (`.energy`) USING HASH\n"
        ") ENGINE=Memory\n"
        "AS\n"
        "(\n"
        "   SELECT\n"
        "      INTERVAL(LOG10(Energy), %100:energyest) AS `.energy`,\n"
        "      COUNT(*) AS CountN,\n"
        "      SUM(    (%101:spectrum)/pow(Energy, SpectralIndex)   ) AS CountW,\n"   // Contents is: CountW
        "      SUM(POW((%101:spectrum)/pow(Energy, SpectralIndex),2)) AS CountW2\n"   // Error    is: SQRT(CountW2)
        "   FROM\n"
        "      MonteCarloFiles\n"
        "   LEFT JOIN\n"
        "      factmc.RunInfoMC USING (FileId)\n"
        "   LEFT JOIN\n"
        "      factmc.OriginalMC USING (FileId)\n"
        "   GROUP BY\n"
        "      `.energy`\n"
        "   ORDER BY\n"
        "      `.energy`\n"
        ")";

    query7.parse();
    for (auto it=env.cbegin(); it!=env.cend(); it++)
        query7.template_defaults[it->first.c_str()] = it->second.c_str();

    //query7.template_defaults["list"]      = listmc.c_str();
    query7.template_defaults["energyest"] = stre.c_str();
    query7.template_defaults["spectrum"]  = spectrum.c_str();

    if (print_queries)
        PrintQuery(query7.str());

    qlog << query7 << ";\n" << endl;
    if (connection.connected())
    {
        cout << query7.execute().info() << endl;
        ShowWarnings(connection);
        Dump(flog, connection, "SimulatedSpectrum");

        const auto sec7 = Time().UnixTime()-start7.UnixTime();
        cout << "Execution time: " << sec7 << "s\n";


        const mysqlpp::StoreQueryResult res7 = connection.query("SELECT * FROM SimulatedSpectrum").store();

        //cout << "Received rows: " << res7.num_rows() << '\n' << endl;

        //       123456789|123456789|123456789|123456789|123456789|123456789|123456789
        cout << "       Bin CountW           CountW2" << endl;
        const auto bins = binning_esim.vec();
        for (auto ir=res7.cbegin(); ir!=res7.cend(); ir++)
        {
            const mysqlpp::Row &row = *ir;

            const uint32_t bin = row[".energy"];
            cout << setw(5) << bins[bin] << ": " << setw(10) << row["CountW"] << " " << setw(10) << row["CountW2"] << '\n';
        }
        cout << endl;
    }


    // -------------------------------------------------------------------
    // --------------------------- 8: Spectrum ---------------------------

    cout << separator("Spectrum") << '\n';

    Time start8;

    /* 02:calculate-spectrum.sql */

    mysqlpp::Query query8(&connection);
    query8 <<
        // SELECT SUM(CountN) AS CountN, SUM(OnTime) AS OnTime FROM ThetaHist
        "CREATE TEMPORARY TABLE Spectrum\n"
        "(\n"
        "   `.energy`      SMALLINT UNSIGNED NOT NULL COMMENT 'Bin Index [Energy]' PRIMARY KEY,\n"
        "   lo             DOUBLE            NOT NULL COMMENT 'Lower edge of energy bin in lg(E/GeV)',\n"
        "   hi             DOUBLE            NOT NULL COMMENT 'Upper edge of energy bin in lg(E/GeV)',\n"
        "   `Signal`       DOUBLE            NOT NULL COMMENT 'Number of signal events',\n"
        "   `Background`   DOUBLE            NOT NULL COMMENT 'Average number of background events',\n"
        "   `Excess`       DOUBLE            NOT NULL COMMENT 'Number of excess events',\n"
        "   ErrSignal      DOUBLE            NOT NULL COMMENT 'Poisson error on number of signal events',\n"
        "   ErrBackground  DOUBLE            NOT NULL COMMENT 'Poisson error on number of background events',\n"
        "   `ErrExcess`    DOUBLE            NOT NULL COMMENT 'Error of excess events',\n"
        "   `Significance` DOUBLE            NOT NULL COMMENT 'Li/Ma sigficance',\n"
        "   `ExcessN`      DOUBLE            NOT NULL COMMENT 'Number of excess events in simulated data',\n"
        "   `ExcessW`      DOUBLE            NOT NULL COMMENT 'Weighted number of excess events in simulated data',\n"
        "   `ErrExcessN`   DOUBLE            NOT NULL COMMENT 'Error or number of excess events in simulated data',\n"
        "   `ErrExcessW`   DOUBLE            NOT NULL COMMENT 'Error of weighted number of excess events in simulated data',\n"
        "   SignalW        DOUBLE            NOT NULL COMMENT 'Weighted number of signal events in simulated data',\n"
        "   BackgroundW    DOUBLE            NOT NULL COMMENT 'Weighted number of background events in simulated data',\n"
        "   ErrSignalW     DOUBLE            NOT NULL COMMENT 'Error of weighted number of signal events in simulated data',\n"
        "   ErrBackgroundW DOUBLE            NOT NULL COMMENT 'Error of weighted number of background events in simulated data',\n"
        "   Flux           DOUBLE            NOT NULL COMMENT 'dN/dA/dt [cm^-2 s-^1 TeV^-1]',\n"
        "   ErrFlux        DOUBLE            NOT NULL COMMENT 'dN/dA/dt [cm^-2 s-^1 TeV^-1]',\n"
        "   Bias           DOUBLE            NOT NULL COMMENT 'Energy Bias, average residual in lg(E)',\n"
        "   Resolution     DOUBLE            NOT NULL COMMENT 'Energy resolution, standard divation of residual in lg(E)',\n"
        "   EfficiencyN    DOUBLE            NOT NULL COMMENT 'Simulated cut efficiency (weighted)',\n"
        "   EfficiencyW    DOUBLE            NOT NULL COMMENT 'Simulated cut efficiency (unweighted)',\n"
        "   ErrEfficiencyN DOUBLE            NOT NULL COMMENT 'Error of simulated cut efficiency (weighted)',\n"
        "   ErrEfficiencyW DOUBLE            NOT NULL COMMENT 'Error of simulated cut efficiency (unweighted)'\n"
        ") ENGINE=Memory\n"
        "AS\n"
        "(\n"
        "   WITH ThetaSums AS\n"
        "   (\n"
        "      SELECT\n"
        "         SUM(CountN) AS CountSim,\n"
        "         SUM(OnTime) AS ObsTime\n"
        "      FROM\n"
        "         ThetaHist\n"
        "   ),\n"
        "   Area AS\n"
        "   (\n"
        "      SELECT\n"
        "         POW(MinImpactHi,2)*PI() AS Area\n"
        "      FROM\n"
        "         MonteCarloArea\n"
        "   ),\n"
        "   ResultMC AS\n" // Group AnalysisMC by EnergyEst Binning
        "   (\n"
        "      SELECT\n"
        "         `.energyest`             AS `.energy`,\n"
        "         ANY_VALUE(SignalW)       AS SignalW,\n"
        "         ANY_VALUE(SignalW2)      AS SignalW2,\n"
        "         ANY_VALUE(BackgroundW)   AS BackgroundW,\n"
        "         ANY_VALUE(BackgroundW2)  AS BackgroundW2,\n"
        "         ANY_VALUE(SignalN)       AS SignalN,\n"
        "         ANY_VALUE(BackgroundN)   AS BackgroundN,\n"
        "         ANY_VALUE(ExcessW)       AS ExcessW,\n"
        "         ANY_VALUE(ExcessN)       AS ExcessN,\n"
        "         ANY_VALUE(ErrExcessW)    AS ErrExcessW,\n"
        "         ANY_VALUE(ErrExcessN)    AS ErrExcessN,\n"
        "         ANY_VALUE(BiasEst)       AS Bias,\n"
        "         ANY_VALUE(ResolutionEst) AS Resolution\n"
        "      FROM\n"
        "         AnalysisMC\n"
        "      GROUP BY\n"
        "         `.energy`\n"
        "      ORDER BY\n"
        "         `.energy`\n"
        "   )\n"
        "   SELECT\n"
        "      `.energy`, lo, hi,\n"  // Scale for Theta-Weights
        "      `Signal`, `Background`/5 AS `Background`, `Excess`, `ErrExcess`, `Significance`,\n"
        "      SQRT(`Signal`)         AS ErrSignal,\n"
        "      SQRT(`SignalW2`)       AS ErrSignalW,\n"
        "      SQRT(`Background`)/5   AS ErrBackground,\n"
        "      SQRT(`BackgroundW2`)/5 AS ErrBackgroundW,\n"
        "      ExcessN, ExcessW, ErrExcessN, ErrExcessW, SignalW, BackgroundW,\n"
        "      AnalysisData.Excess/ResultMC.ExcessW*SimulatedSpectrum.CountW * 1000/(POW(10,hi)-POW(10,lo)) /Area/ObsTime / CountSim*ObsTime AS Flux,\n"
        "      AnalysisData.Excess/ResultMC.ExcessW*SimulatedSpectrum.CountW * 1000/(POW(10,hi)-POW(10,lo)) /Area/ObsTime / CountSim*ObsTime\n"
        "         * SQRT(\n"
        "             + POW(AnalysisData.ErrExcess / AnalysisData.Excess, 2)\n"
        "             + POW(ResultMC.ErrExcessW    / ResultMC.ExcessW,    2)\n"
        "             + SimulatedSpectrum.CountW2  / POW(SimulatedSpectrum.CountW,2)\n"
        "           ) AS ErrFlux,\n"
        "      Bias,\n"
        "      Resolution,\n"
        "      ResultMC.ExcessW/SimulatedSpectrum.CountW * CountSim/ObsTime AS EfficiencyW,\n"
        "      ResultMC.ExcessN/SimulatedSpectrum.CountN AS EfficiencyN,\n"
        "      ( POW(ResultMC.ErrExcessW/ResultMC.ExcessW, 2) + POW(SQRT(SimulatedSpectrum.CountW2)/SimulatedSpectrum.CountW, 2) )\n"
        "         * POW(ResultMC.ExcessW/SimulatedSpectrum.CountW, 2) * CountSim/ObsTime AS ErrEfficiencyW,\n"
        "      ( POW(ResultMC.ErrExcessN, 2) + POW(ResultMC.ExcessN, 2)/SimulatedSpectrum.CountN)/POW(SimulatedSpectrum.CountN, 2) AS ErrEfficiencyN\n"
        "   FROM\n"
        "      AnalysisData\n"
        "   INNER JOIN\n"
        "      ResultMC USING(`.energy`)\n"
        "   INNER JOIN\n"
        "      SimulatedSpectrum USING(`.energy`)\n"
        "   INNER JOIN\n"
        "      BinningEnergyEst ON `.energy`=bin\n"
        "   CROSS JOIN\n"
        "      ThetaSums, Area\n"
        "   WHERE\n"
        "      AnalysisData.Excess>0\n"
        "   ORDER BY\n"
        "      `.energy`\n"
        ")"
        ;

    // [ Sa^2/a^2 + Sb^2/b^2 ] * a^2/b^2
    // [ (sc^2)/c^2+(sb^2)/b^2+(sa^2)/a^2 ]  * a^2*b^2/c^2

    query8.parse();
    for (auto it=env.cbegin(); it!=env.cend(); it++)
        query8.template_defaults[it->first.c_str()] = it->second.c_str();

    //query8.template_defaults["area"] = area;
    //query8.template_defaults["ontime"] = resX[0]["OnTime"].data();
    //query8.template_defaults["count"]  = resX[0]["CountN"].data();

    if (print_queries)
        PrintQuery(query8.str());

    qlog << query8 << ";\n" << endl;
    if (connection.connected())
    {
        cout << query8.execute().info() << endl;
        ShowWarnings(connection);
        Dump(flog, connection, "Spectrum");

        const auto sec8 = Time().UnixTime()-start8.UnixTime();
        cout << "Execution time: " << sec8 << "s\n";


        const mysqlpp::StoreQueryResult res8 = connection.query("SELECT * FROM Spectrum").store();

        cout << "  Bin  Flux                   Error" << endl;
        const auto bins = binning_eest.vec();
        for (auto ir=res8.cbegin(); ir!=res8.cend(); ir++)
        {
            const mysqlpp::Row &row = *ir;

            const uint32_t bin = row[".energy"];
            cout << setw(5) << bins[bin] << ": " << setw(10) << row["Flux"] << " " << setw(10) << row["ErrFlux"] << '\n';
        }
        cout << endl;

        // --------------------------------------------------------------------------

        // Crab Nebula: 1 TeV: 3e-7  /  m^2 / s / TeV
        // Crab Nebula: 1 TeV: 3e-11 / cm^2 / s / TeV

        sindent indentm(mlog);

        mlog << "void spectrum()\n";
        mlog << "{\n" << indent(4);
        mlog << "// Energy Spectrum (e.g. Crab: 3e-11 [cm^-2 s-^1 TeV^-1] @ 1TeV)\n";
        mlog << "TGraphErrors g;\n";
        mlog << "g.SetNameTitle(\"Spectrum\", \"Energy Spectrum\");\n";
        for (auto ir=res8.cbegin(); ir!=res8.cend(); ir++)
        {
            // This is not efficient but easier to read and safer
            const mysqlpp::Row &row = *ir;

            const double hi     = row["hi"];
            const double lo     = row["lo"];
            const double center = (hi+lo)/2;
            const double flux   = row["Flux"];
            const double error  = row["ErrFlux"];

            mlog << "g.SetPoint(g.GetN(), " << setw(8) << center << ", " << flux << ");\n";
            mlog << "g.SetPointError(g.GetN()-1, 0, " << error << ");\n";
        }
        mlog << "g.GetXaxis()->SetTitle(\"lg(E/TeV)\");\n";
        mlog << "g.GetYaxis()->SetTitle(\"dN/dE [cm^{-2} s^{-1} TeV^{-1}]\");";
        mlog << endl;
    }

    WriteHistogram(connection, {
             .name     = "Signal",
             .title    = "Signal",
             .dir      = "Eest/Measurement",
             .binningx = binning_eest,
             .table    = "Spectrum",
             .x        = ".energy",
             .v        = "Signal",
             .err      = "ErrSignal",
             .axisx    = "Energy lg(E_{est}/GeV)",
             .axisy    = "Counts"
         });

    WriteHistogram(connection, {
             .name     = "Background",
             .title    = "Background",
             .dir      = "Eest/Measurement",
             .binningx = binning_eest,
             .table    = "Spectrum",
             .x        = ".energy",
             .v        = "Background",
             .err      = "ErrBackground",
             .axisx    = "Energy lg(E_{est}/GeV)",
             .axisy    = "Counts"
         });

    WriteHistogram(connection, {
             .name     = "Excess",
             .title    = "Excess",
             .dir      = "Eest/Measurement",
             .binningx = binning_eest,
             .table    = "Spectrum",
             .x        = ".energy",
             .v        = "Excess",
             .err      = "ErrExcess",
             .axisx    = "Energy lg(E_{est}/GeV)",
             .axisy    = "Signal - Background (Counts)"
         });

    WriteHistogram(connection, {
             .name     = "SignalW",
             .title    = "SignalW",
             .dir      = "Eest/Simulation/Weighted",
             .binningx = binning_eest,
             .table    = "Spectrum",
             .x        = ".energy",
             .v        = "SignalW",
             .err      = "ErrSignalW",
             .axisx    = "Energy lg(E_{est}/GeV)",
             .axisy    = "Weighted"
         });

    WriteHistogram(connection, {
             .name     = "BackgroundW",
             .title    = "BackgroundW",
             .dir      = "Eest/Simulation/Weighted",
             .binningx = binning_eest,
             .table    = "Spectrum",
             .x        = ".energy",
             .v        = "BackgroundW",
             .err      = "ErrBackgroundW",
             .axisx    = "Energy lg(E_{est}/GeV)",
             .axisy    = "Weighted"
         });

    WriteHistogram(connection, {
             .name     = "ExcessW",
             .title    = "ExcessW",
             .dir      = "Eest/Simulation/Weighted",
             .binningx = binning_eest,
             .table    = "Spectrum",
             .x        = ".energy",
             .v        = "ExcessW",
             .err      = "ErrExcessW",
             .axisx    = "Energy lg(E_{est}/GeV)",
             .axisy    = "Signal - Background (Weighted)"
         });

    WriteHistogram(connection, {
             .name     = "Significance",
             .title    = "Significance",
             .dir      = "Eest/Measurement",
             .binningx = binning_eest,
             .table    = "Spectrum",
             .x        = ".energy",
             .v        = "Significance",
             .axisx    = "Energy lg(E_{est}/GeV)",
             .axisy    = "Li/Ma Significance"
         });

    WriteHistogram(connection, {
             .dir      = "Eest",
             .name     = "Bias",
             .title    = "Energy Bias",
             .binningx = binning_eest,
             .table    = "Spectrum",
             .x        = ".energy",
             .v        = "Bias",
             .err      = "Resolution",
             .axisx    = "Energy lg(E_{sim}/GeV)",
             .axisy    = "lg(E_{est}/E_{sim})",
             .stats    = false
         });

    WriteHistogram(connection, {
             .dir      = "Eest",
             .name     = "Resolution",
             .title    = "Energy Resolution",
             .binningx = binning_eest,
             .table    = "Spectrum",
             .x        = ".energy",
             .v        = "Resolution",
             .axisx    = "Energy lg(E_{sim}/GeV)",
             .axisy    = "#sigma(lg(E_{est}/E_{sim}))",
             .stats    = false
         });

    WriteHistogram(connection, {
             .dir      = "Eest/Efficiency",
             .name     = "EfficiencyN",
             .title    = "Efficiency (Counts)",
             .binningx = binning_eest,
             .table    = "Spectrum",
             .x        = ".energy",
             .v        = "EfficiencyN",
             .err      = "ErrEfficiencyN",
             .axisx    = "Energy lg(E_{sim}/GeV)",
             .axisy    = "Efficiency"
         });

    WriteHistogram(connection, {
             .dir      = "Eest/Efficiency",
             .name     = "EfficiencyW",
             .title    = "Efficiency (Weighted)",
             .binningx = binning_eest,
             .table    = "Spectrum",
             .x        = ".energy",
             .v        = "EfficiencyW",
             .err      = "ErrEfficiencyW",
             .axisx    = "Energy lg(E_{sim}/GeV)",
             .axisy    = "Efficiency"
         });

    WriteHistogram(connection, {
             .name     = "Spectrum",
             .title    = "Differential Energy Spectrum",
             .binningx = binning_eest,
             .table    = "Spectrum",
             .x        = ".energy",
             .v        = "Flux",
             .err      = "ErrFlux",
             .axisx    = "Energy lg(E/GeV)",
             .axisy    = "dN/dE [cm^{-2} s^{-1} TeV^{-1}]",
             .stats    = false
         });


    // -------------------------------------------------------------------
    // -------------------------- 9: Threshold ---------------------------

    cout << separator("Threshold") << '\n';

    Time start9;

    /* 02:calculate-threshold.sql */

    // This query gets the analysis results versus Simulated Energy from the combined table
    mysqlpp::Query query9(&connection);
    query9 <<
        // SELECT SUM(CountN) AS CountN, SUM(OnTime) AS OnTime FROM ThetaHist
        "CREATE TEMPORARY TABLE Threshold ENGINE=Memory AS\n"
        "(\n"
        "   WITH ThetaSums AS\n"
        "   (\n"
        "      SELECT\n"
        "         SUM(CountN) AS CountSim,\n"
        "         SUM(OnTime) AS ObsTime\n"
        "      FROM\n"
        "         ThetaHist\n"
        "   ),\n"
        "   Area AS\n"
        "   (\n"
        "      SELECT\n"
        "         POW(MinImpactHi,2)*PI() AS Area\n"
        "      FROM\n"
        "         MonteCarloArea\n"
        "   ),\n"
        "   ResultMC AS\n" // Group AnalysisMC by EnergySim Binning
        "   (\n"
        "      SELECT\n"
        "         `.energysim`             AS `.energy`,\n"
        "         ANY_VALUE(ThresholdW)    AS ThresholdW,\n"
        "         ANY_VALUE(ThresholdW2)   AS ThresholdW2,\n"
        "         ANY_VALUE(ThresholdN)    AS ThresholdN,\n"
        "         ANY_VALUE(BiasSim)       AS Bias,\n"
        "         ANY_VALUE(ResolutionSim) AS Resolution\n"
        "      FROM\n"
        "         AnalysisMC\n"
        "      GROUP BY\n"
        "         `.energy`\n"
        "   )\n"
        "   SELECT\n"
        "      `.energy`, lo, hi,\n"
        "      ThresholdW,\n"
        "      SQRT(ThresholdW2) AS ErrThresholdW,\n"
        "      ThresholdN,\n"
        "      SQRT(ThresholdN)  AS ErrThresholdN,\n"         // Scale for Theta-Weights
        "      ThresholdW        * 1000/(POW(10,hi)-POW(10,lo)) / Area / CountSim*ObsTime AS Flux,\n"
        "      SQRT(ThresholdW2) * 1000/(POW(10,hi)-POW(10,lo)) / Area / CountSim*ObsTime AS ErrFlux,\n"
        "      Bias,\n"
        "      Resolution\n"
        "   FROM\n"
        "      ResultMC\n"
        "   INNER JOIN\n"
        "      BinningEnergySim ON `.energy`=bin\n"
        "   CROSS JOIN\n"
        "      ThetaSums, Area\n"
        "   WHERE\n"
        "      ThresholdW>0 AND ThresholdW2>0\n"
        "   ORDER BY\n"
        "      `.energy`\n"
        ")";

    query9.parse();
    for (auto it=env.cbegin(); it!=env.cend(); it++)
        query9.template_defaults[it->first.c_str()] = it->second.c_str();

    //query9.template_defaults["area"] = area;
    //query9.template_defaults["ontime"] = resX[0]["OnTime"].data();
    //query9.template_defaults["count"]  = resX[0]["CountN"].data();

    if (print_queries)
        PrintQuery(query9.str());

    qlog << query9 << ";\n" << endl;
    if (connection.connected())
    {
        cout << query9.execute().info() << endl;
        ShowWarnings(connection);
        Dump(flog, connection, "Threshold");

        const auto sec9 = Time().UnixTime()-start9.UnixTime();
        cout << "Execution time: " << sec9 << "s\n";

        // --------------------------------------------------------------------------

        const mysqlpp::StoreQueryResult res9 = connection.query("SELECT * FROM Threshold").store();

        sindent indentm(mlog, 4);

        mlog << '\n';
        mlog << "// Energy Threshold\n";
        mlog << "TGraphErrors g;\n";
        mlog << "g.SetNameTitle(\"Threshold\", \"Energy Threshold\");\n";
        for (auto ir=res9.cbegin(); ir!=res9.cend(); ir++)
        {
            // This is not efficient but easier to read and safer
            const mysqlpp::Row &row = *ir;

            const double hi     = row["hi"];
            const double lo     = row["lo"];
            const double center = (hi+lo)/2;
            const double width  = pow(10, hi)-pow(10, lo);
            const double flux   = row["Flux"]   /width;
            const double error  = row["ErrFlux"]/width;

            mlog << "g.SetPoint(g.GetN(), " << setw(8) << center << ", " << flux << ");\n";
            mlog << "g.SetPointError(g.GetN()-1, 0, " << error << ");\n";
        }
        mlog << "g.GetXaxis()->SetTitle(\"lg(E/TeV)\");\n";
        mlog << "g.GetYaxis()->SetTitle(\"dN/dE [cm^{-2} s^{-1} TeV^{-1}]\");\n";
        mlog << indent(0) << "}" << endl;
    }

    WriteHistogram(connection, {
             .name     = "SimExcessW",
             .title    = "Weighted Simulated Excess",
             .dir      = "Esim/Simulation/Weighted",
             .binningx = binning_esim,
             .table    = "Threshold",
             .x        = ".energy",
             .v        = "ThresholdW",
             .err      = "ErrThresholdW",
             .axisx    = "Energy lg(E_{sim}/GeV)",
             .axisy    = "Weighted Counts"
         });

    WriteHistogram(connection, {
             .name     = "SimExcessN",
             .title    = "Simulated Excess",
             .dir      = "Esim/Simulation/Counts",
             .binningx = binning_esim,
             .table    = "Threshold",
             .x        = ".energy",
             .v        = "ThresholdN",
             .err      = "ErrThresholdN",
             .axisx    = "Energy lg(E_{sim}/GeV)",
             .axisy    = "Counts"
         });

    WriteHistogram(connection, {
             .name     = "SimSpectrumW",
             .title    = "Weighted Simulated Excess Spectrum",
             .binningx = binning_esim,
             .dir      = "Esim/Simulation/Weighted",
             .table    = "Threshold",
             .x        = ".energy",
             .v        = "Flux",
             .err      = "ErrFlux",
             .axisx    = "Energy lg(E_{sim}/GeV)",
             .axisy    = "dN/dE [cm^{-2} s^{-1} TeV^{-1}]"
         });

    WriteHistogram(connection, {
             .dir      = "Esim",
             .name     = "Bias",
             .title    = "Energy Bias",
             .binningx = binning_esim,
             .table    = "Threshold",
             .x        = ".energy",
             .v        = "Bias",
             .err      = "Resolution",
             .axisx    = "Energy lg(E_{sim}/GeV)",
             .axisy    = "lg(E_{est}/E_{sim})",
             .stats    = false
         });

    WriteHistogram(connection, {
             .dir      = "Esim",
             .name     = "Resolution",
             .title    = "Energy Resolution",
             .binningx = binning_esim,
             .table    = "Threshold",
             .x        = ".energy",
             .v        = "Resolution",
             .axisx    = "Energy lg(E_{sim}/GeV)",
             .axisy    = "#sigma(lg(E_{est}/E_{sim}))",
             .stats    = false
         });


    // -------------------------------------------------------------------
    // -------------------------- 10: Migration --------------------------

    cout << separator("Migration") << '\n';

    Time start10;

    /* 02:obtain-migration-matrix.sql */

    // This query gets the analysis results versus Estimated Energy from the combined table
    mysqlpp::Query query10(&connection);
    query10 <<
        "CREATE TEMPORARY TABLE Migration ENGINE=Memory AS\n"
        "(\n"
        "   SELECT\n"
        "      `.energyest`,\n"
        "      `.energysim`,\n"
        "      BinningEnergySim.lo   AS EsimLo,\n"
        "      BinningEnergySim.hi   AS EsimHi,\n"
        "      BinningEnergyEst.lo   AS EestLo,\n"
        "      BinningEnergyEst.hi   AS EestHi,\n"
        "      ANY_VALUE(MigrationW) AS MigrationW,\n"
        "      ANY_VALUE(MigrationN) AS MigrationN\n"
        //     FIXME: Errors
        "   FROM\n"
        "      AnalysisMC\n"
        "   INNER JOIN\n"
        "      BinningEnergyEst ON `.energyest`=BinningEnergyEst.bin\n"
        "   INNER JOIN\n"
        "      BinningEnergySim ON `.energysim`=BinningEnergySim.bin\n"
        "   GROUP BY\n"
        "      `.energyest`, `.energysim`\n"
        "   ORDER BY\n"
        "      `.energyest`, `.energysim`\n"
        ")";

    if (print_queries)
        PrintQuery(query10.str());

    qlog << query10 << ";\n" << endl;
    if (connection.connected())
    {
        cout << query10.execute().info() << endl;
        ShowWarnings(connection);
        Dump(flog, connection, "Migration");

        const auto sec10 = Time().UnixTime()-start10.UnixTime();
        cout << "Execution time: " << sec10 << "s\n\n";
    }

    WriteHistogram(connection, {
             .name     = "MigrationN",
             .title    = "Energy Migration",
             .binningx = binning_esim,
             .binningy = binning_eest,
             .table    = "Migration",
             .x        = ".energysim",
             .y        = ".energyest",
             .v        = "MigrationN",
             .axisx    = "Energy lg(E_{sim}/GeV)",
             .axisy    = "Energy lg(E_{est}/GeV)",
             .axisz    = "Counts",
             .stats    = false
         });

    WriteHistogram(connection, {
             .name     = "MigrationW",
             .title    = "Energy Migration",
             .binningx = binning_esim,
             .binningy = binning_eest,
             .table    = "Migration",
             .x        = ".energysim",
             .y        = ".energyest",
             .v        = "MigrationW",
             .axisx    = "Energy lg(E_{sim}/GeV)",
             .axisy    = "Energy lg(E_{est}/GeV)",
             .axisz    = "Weighted Counts",
             .stats    = false
         });



    // -------------------------------------------------------------------
    // --------------------------- 11: Summary ---------------------------

    cout << separator("Summary") << '\n';
    const auto sec = Time().UnixTime()-start.UnixTime();
    cout << "Total execution time: " << sec << "s\n" << endl;

    return 0;
}