Changeset 19895


Ignore:
Timestamp:
Dec 14, 2019, 7:39:38 PM (4 months ago)
Author:
tbretz
Message:
Completele redone the logic, now the program produces all kind of binning and combinations that one can think of (at least almost) and also automatically calculates upper limits. Still missing: integral spectra.
File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/FACT++/src/spectrum.cc

    r19889 r19895  
    1717#include "TH1.h"
    1818#include "TH2.h"
     19#include "TRolke.h"
     20#include "TFeldmanCousins.h"
    1921#endif
    2022
     
    130132#endif
    131133         , "Database link as in\n\tuser:password@server[:port]/database[?compress=0|1].")
    132         ("out,o",            var<string>(conf.GetName()), "Defines the prefix (with path) of the output files.")
    133         ("dry-run",          po_switch(),   "Only write the queries to the query.log file. Internally overwrites uri with an empty string.")
     134        ("out,o", var<string>(conf.GetName()), "Defines the prefix (with path) of the output files.")
     135        ("confidence-level,c", var<double>(0.99), "Confidence level for the calculation of the upper limits.")
     136        ("feldman-cousins", po_bool(), "Calculate Feldman-Cousins ULs (slow and only minor difference to Rolke).")
    134137        ;
    135138
    136139    po::options_description binnings("Binnings");
    137140    binnings.add_options()
    138         ("theta",     var<Binning>(Binning(90, 0, 90)), "Add equidistant bins in theta (degrees). Syntax: N,lo,hi (Number N of equidistant bins between lo and hi)")
    139         ("theta-bin", vars<double>(),                   "Add a bin-edge to the theta binning (degree)")
    140         ("esim",      var<Binning>(Binning(15, 2, 5)),  "Add equidistant bins in log10 simulated energy. Syntax: N,lo,hi (Number N of equidistant bins between lo and hi)")
    141         ("esim-bin",  vars<double>(),                   "Add a bin-edge to the binnnig in log10 simulated enegry")
    142         ("eest",      var<Binning>(Binning(15, 2, 5)),  "Add equidistant bins in log10 estimated energy. Syntax: N,lo,hi (Number N of equidistant bins between lo and hi)")
    143         ("eest-bin", vars<double>(),                   "Add a bin-edge to the binning in log10 estimated enegry")
     141        ("theta",             var<Binning>(Binning(90, 0, 90)), "Add equidistant bins in theta (degrees). Syntax: N,lo,hi (Number N of equidistant bins between lo and hi)")
     142        ("theta-bin",         vars<double>(),                   "Add a bin-edge to the theta binning (degree)")
     143        ("energy-dense",      var<Binning>(Binning(30, 2, 5)),  "Add equidistant bins in log10 simulated energy. Syntax: N,lo,hi (Number N of equidistant bins between lo and hi)")
     144        ("energy-dense-bin",  vars<double>(),                   "Add a bin-edge to the binnnig in log10 simulated enegry")
     145        ("energy-sparse",     var<Binning>(Binning(15, 2, 5)),  "Add equidistant bins in log10 estimated energy. Syntax: N,lo,hi (Number N of equidistant bins between lo and hi)")
     146        ("energy-sparse-bin", vars<double>(),                   "Add a bin-edge to the binning in log10 estimated enegry")
    144147        ;
    145148
     
    147150    queries.add_options()
    148151        ("analysis",    var<string>("analysis.sql"),   "File with the analysis query. A default file is created automatically in the <prefix> directory it does not exist.")
    149         ("data",        var<string>("data.sql"),       "File with the query which creates the data summary. A default file is created automatically in the <prefix> directory it does not exist.")
    150         ("simulation",  var<string>("simulation.sql"), "File with the query which creates the simulation summary. A default file is created automatically in the <prefix> directory it does not exist.")
     152        //("data",        var<string>("data.sql"),       "File with the query which creates the data summary. A default file is created automatically in the <prefix> directory it does not exist.")
     153        //("simulation",  var<string>("simulation.sql"), "File with the query which creates the simulation summary. A default file is created automatically in the <prefix> directory it does not exist.")
    151154        ;
    152155
     
    162165    po::options_description debug("Debug options");
    163166    debug.add_options()
    164         ("print-connection", po_switch(),       "Print database connection information")
     167        ("dry-run",          po_bool(),         "Only write the queries to the query.log file. Internally overwrites uri with an empty string.")
     168        ("print-connection", po_bool(),         "Print database connection information")
    165169#ifdef HAVE_HIGHLIGHT
    166         ("print-queries",    po_switch(),       "Print all queries to the console. They are automatically highlighted.")
     170        ("print-queries",    po_bool(),         "Print all queries to the console. They are automatically highlighted.")
    167171#else
    168         ("print-queries",    po_switch(),       "Print all queries to the console. (For highlighting recompile with 'libsource-highlight-dev' installed)")
     172        ("print-queries",    po_bool(),         "Print all queries to the console. (For highlighting recompile with 'libsource-highlight-dev' installed)")
    169173#endif
     174        ("mc-only",          po_bool(),         "Do not run an data realated queries (except observation times)")
    170175        ("verbose,v",        var<uint16_t>(1),  "Verbosity (0: quiet, 1: default, 2: more, 3, ...)")
    171176        ;
     
    327332    srchilite::SourceHighlight sourceHighlight("esc256.outlang");
    328333    sourceHighlight.setStyleFile("esc256.style");
     334    sourceHighlight.setGenerateLineNumbers();
     335    sourceHighlight.setLineNumberDigits(3);
     336    //sourceHighlight.setLineNumberPad(' ')
    329337    sourceHighlight.highlight(qstr, cout, "sql.lang");
    330338    cout << endl;
     
    341349        "(\n"
    342350        "   bin SMALLINT UNSIGNED NOT NULL,\n"
    343         "   lo FLOAT NOT NULL,\n"
    344         "   hi FLOAT NOT NULL,\n"
     351        "   lo DOUBLE NOT NULL,\n"
     352        "   hi DOUBLE NOT NULL,\n"
    345353        "   PRIMARY KEY (bin) USING HASH\n"
    346354        ")";
     
    349357    if (connection.connected())
    350358        query0.execute();
    351 
    352     //connection.query("ALTER TABLE Binning"+name+" AUTO_INCREMENT=0").execute();
    353359
    354360    mysqlpp::Query query1(&connection);
     
    426432};
    427433
     434#ifdef HAVE_ROOT
     435
     436TH1 *CreateHistogram(const Histogram &hist)
     437{
     438    const char *name = hist.name.empty() ? hist.v.c_str() : hist.name.c_str();
     439
     440    cout << "Creating Histogram '" << hist.dir << "/" << name << "'" << endl;
     441
     442    const auto vecx = hist.binningx.vec();
     443    const auto vecy = hist.binningy.vec();
     444
     445    TH1 *h = 0;
     446
     447    if (hist.y.empty())
     448    {
     449        h = hist.binningx.equidist ?
     450            new TH1D(name, hist.title.c_str(), vecx.size()-1, vecx.front(), vecx.back()) :
     451            new TH1D(name, hist.title.c_str(), vecx.size()-1, vecx.data());
     452    }
     453    else
     454    {
     455        if (hist.binningy.equidist)
     456        {
     457            h = hist.binningx.equidist ?
     458                new TH2D(name, hist.title.c_str(), vecx.size()-1, vecx.front(), vecx.back(), vecy.size()-1, vecy.front(), vecy.back()) :
     459                new TH2D(name, hist.title.c_str(), vecx.size()-1, vecx.data(), vecy.size()-1, vecy.front(), vecy.back());
     460        }
     461        else
     462        {
     463            h = hist.binningx.equidist ?
     464                new TH2D(name, hist.title.c_str(), vecx.size()-1, vecx.front(), vecx.back(), vecy.size()-1, vecy.front(), vecy.back()) :
     465                new TH2D(name, hist.title.c_str(), vecx.size()-1, vecx.data(), vecy.size()-1, vecy.front(), vecy.back());
     466        }
     467    }
     468
     469    h->SetXTitle(hist.axisx.c_str());
     470    h->SetYTitle(hist.axisy.c_str());
     471    h->SetZTitle(hist.axisz.c_str());
     472
     473    h->SetMarkerStyle(kFullDotMedium);
     474
     475    h->SetStats(hist.stats);
     476
     477    return h;
     478}
     479
     480void WriteHistogram(TH1 *h, const string &directory)
     481{
     482    gFile->cd();
     483    TDirectory *dir = gFile->GetDirectory(directory.c_str());
     484    if (dir)
     485        dir->cd();
     486    else
     487    {
     488        gFile->mkdir(directory.c_str());
     489        gFile->cd(directory.c_str());
     490    }
     491    h->Write();
     492}
     493#endif
     494
    428495void WriteHistogram(Database &connection, const Histogram &hist)
    429496{
     
    431498    if (!connection.connected())
    432499        return;
     500
     501    TH1 *h = CreateHistogram(hist);
    433502
    434503    mysqlpp::Query query(&connection);
     
    449518    const mysqlpp::StoreQueryResult res = query.store();
    450519
    451     const auto vecx = hist.binningx.vec();
    452     const auto vecy = hist.binningy.vec();
    453 
    454     TH1 *h = 0;
    455 
    456     if (hist.y.empty())
    457     {
    458         h = hist.binningx.equidist ?
    459             new TH1D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.front(), vecx.back()) :
    460             new TH1D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.data());
    461     }
    462     else
    463     {
    464         if (hist.binningy.equidist)
    465         {
    466             h = hist.binningx.equidist ?
    467                 new TH2D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.front(), vecx.back(), vecy.size()-1, vecy.front(), vecy.back()) :
    468                 new TH2D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.data(), vecy.size()-1, vecy.front(), vecy.back());
    469         }
    470         else
    471         {
    472             h = hist.binningx.equidist ?
    473                 new TH2D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.front(), vecx.back(), vecy.size()-1, vecy.front(), vecy.back()) :
    474                 new TH2D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.data(), vecy.size()-1, vecy.front(), vecy.back());
    475         }
    476     }
    477 
    478     h->SetXTitle(hist.axisx.c_str());
    479     h->SetYTitle(hist.axisy.c_str());
    480     h->SetZTitle(hist.axisz.c_str());
    481 
    482     h->SetMarkerStyle(kFullDotMedium);
    483 
    484     h->SetStats(hist.stats);
    485 
    486520    for (auto ir=res.cbegin(); ir!=res.cend(); ir++)
    487521    {
    488522        const auto &row = *ir;
    489 
    490         const uint32_t x = row["X"];
    491         const double   v = row["V"];
    492         if (x>uint32_t(h->GetNbinsX()+1))
    493             continue;
    494523
    495524        try
    496525        {
    497             const uint32_t y = row["Y"];
    498             if (y>uint32_t(h->GetNbinsY()+1))
     526            const uint32_t x = row["X"];
     527            const double   v = row["V"];
     528            if (x>uint32_t(h->GetNbinsX()+1))
    499529                continue;
    500530
    501             h->SetBinContent(x, y, v);
    502 
    503         }
    504         catch (const mysqlpp::BadFieldName &)
    505         {
    506             h->SetBinContent(x, v);
    507531            try
    508532            {
    509                 h->SetBinError(x, double(row["E"]));
     533                const uint32_t y = row["Y"];
     534                if (y>uint32_t(h->GetNbinsY()+1))
     535                    continue;
     536
     537                h->SetBinContent(x, y, v);
     538
    510539            }
    511540            catch (const mysqlpp::BadFieldName &)
    512541            {
     542                h->SetBinContent(x, v);
     543                try
     544                {
     545                    h->SetBinError(x, double(row["E"]));
     546                }
     547                catch (const mysqlpp::BadFieldName &)
     548                {
     549                }
    513550            }
    514551        }
    515     }
    516 
    517     gFile->cd();
    518     TDirectory *dir = gFile->GetDirectory(hist.dir.c_str());
    519     if (dir)
    520         dir->cd();
    521     else
    522     {
    523         gFile->mkdir(hist.dir.c_str());
    524         gFile->cd(hist.dir.c_str());
    525     }
    526     h->Write();
     552        catch (const mysqlpp::BadConversion &b)
     553        {
     554            cerr << b.what() << endl;
     555        }
     556    }
     557
     558    WriteHistogram(h, hist.dir);
     559    delete h;
     560#endif
     561}
     562
     563void WriteHistogram(Database &connection, const Histogram &hist, const map<size_t, double> &data)
     564{
     565#ifdef HAVE_ROOT
     566    TH1 *h = CreateHistogram(hist);
     567
     568    for (auto ir=data.cbegin(); ir!=data.cend(); ir++)
     569        h->SetBinContent(ir->first, ir->second);
     570
     571    WriteHistogram(h, hist.dir);
    527572    delete h;
    528573#endif
     
    581626    // ----------------------------- Evaluate options --------------------------
    582627
    583     const string   uri     = conf.Get<bool>("dry-run") ? "" : conf.Get<string>("uri");
    584     const string   out     = conf.Get<string>("out");
    585     const uint16_t verbose = conf.Get<uint16_t>("verbose");
     628    const string   uri        = conf.Get<bool>("dry-run") ? "" : conf.Get<string>("uri");
     629    const string   out        = conf.Get<string>("out");
     630    const uint16_t verbose    = conf.Get<uint16_t>("verbose");
     631    const double   confidence = conf.Get<double>("confidence-level");
     632    const bool     feldman    = conf.Get<bool>("feldman-cousins");
    586633
    587634    const bool print_connection = conf.Get<bool>("print-connection");
    588635    const bool print_queries    = conf.Get<bool>("print-queries");
    589 
    590     Binning binning_theta = conf.Get<Binning>("theta")+conf.Vec<double>("theta-bin");
    591     Binning binning_esim  = conf.Get<Binning>("esim");
    592     Binning binning_eest  = conf.Get<Binning>("eest");
     636    const bool mc_only          = conf.Get<bool>("mc-only");
     637
     638    Binning binning_theta  = conf.Get<Binning>("theta")+conf.Vec<double>("theta-bin");
     639    Binning binning_dense  = conf.Get<Binning>("energy-dense");
     640    Binning binning_sparse = conf.Get<Binning>("energy-sparse");
    593641
    594642    cout << '\n';
    595     cout << "Binning 'theta': " << binning_theta.str() << endl;
    596     cout << "Binning 'Esim':  " << binning_esim.str()  << endl;
    597     cout << "Binning 'Eest':  " << binning_eest.str()  << endl;
     643    cout << "Binning 'theta':  " << binning_theta.str() << endl;
     644    cout << "Binning 'dense':  " << binning_dense.str()  << endl;
     645    cout << "Binning 'sparse': " << binning_sparse.str()  << endl;
    598646
    599647    const uint16_t source_key = conf.Get<uint16_t>("source-key");
     
    604652
    605653    cout << "\n";
    606     const string analysis_sql   = CreateResource(conf, "analysis",   "analysis.sql",   RESOURCE(std::string, spectrum_analysis_sql));
    607     const string data_sql       = CreateResource(conf, "data",       "data.sql",       RESOURCE(std::string, spectrum_data_sql));
    608     const string simulation_sql = CreateResource(conf, "simulation", "simulation.sql", RESOURCE(std::string, spectrum_simulation_sql));
     654    const string analysis_sql    = CreateResource(conf, "analysis", "analysis.sql", RESOURCE(std::string, spectrum_analysis_sql));
     655    const string data_sql        = RESOURCE(std::string, spectrum_data_sql);
     656    const string simulation_sql  = RESOURCE(std::string, spectrum_simulation_sql);
     657    const string spectrum_sql    = RESOURCE(std::string, spectrum_spectrum_sql);
     658    const string summary_sim_sql = RESOURCE(std::string, spectrum_summary_sim_sql);
     659    const string summary_est_sql = RESOURCE(std::string, spectrum_summary_est_sql);
    609660    cout << endl;
    610661
    611     const string strzd = binning_theta.list();
    612     const string stre  = binning_eest.list();
    613     const string strth = binning_esim.list();
     662    const string str_theta = binning_theta.list();
     663    const string str_dense  = binning_dense.list();
     664    const string str_sparse = binning_sparse.list();
    614665
    615666    // -------------------------------------------------------------------------
     
    696747    cout << separator("Binnings") << '\n';
    697748
    698     CreateBinning(connection, qlog, binning_theta, "Theta");
    699     CreateBinning(connection, qlog, binning_eest,  "EnergyEst");
    700     CreateBinning(connection, qlog, binning_esim,  "EnergySim");
     749    CreateBinning(connection, qlog, binning_theta,  "Theta");
     750    CreateBinning(connection, qlog, binning_dense,  "Energy_dense");
     751    CreateBinning(connection, qlog, binning_sparse, "Energy_sparse");
    701752
    702753    Dump(flog, connection, "BinningTheta");
    703     Dump(flog, connection, "BinningEnergyEst");
    704     Dump(flog, connection, "BinningEnergySim");
     754    Dump(flog, connection, "BinningEnergy_dense");
     755    Dump(flog, connection, "BinningEnergy_sparse");
    705756
    706757    // -------------------------------------------------------------------
     
    719770        "   FileId INT UNSIGNED NOT NULL,\n"
    720771        "   PRIMARY KEY (FileId) USING HASH\n"
    721         ") ENGINE=Memory\n"
     772        ") ENGINE=MEMORY\n"
    722773        "AS\n"
    723774        "(\n"
     
    770821        "(\n"
    771822        "   `.theta` SMALLINT UNSIGNED NOT NULL,\n"
    772         "   OnTime FLOAT NOT NULL,\n"
     823        "   OnTime DOUBLE NOT NULL,\n"
    773824        "   PRIMARY KEY (`.theta`) USING HASH\n"
    774         ") ENGINE=Memory\n"
     825        ") ENGINE=MEMORY\n"
    775826        "AS\n"
    776827        "(\n"
     
    792843        query2.template_defaults[it->first.c_str()] = it->second.c_str();
    793844
    794     query2.template_defaults["bins"] = strzd.c_str();
     845    query2.template_defaults["bins"] = str_theta.c_str();
    795846
    796847    if (print_queries)
     
    808859    }
    809860
    810     /*
    811 
    812     SELECT
    813        min(first_number)   as first_number,
    814        max(last_number)    as last_number
    815     from
    816     (
    817        select
    818           first_number,
    819           last_number
    820           sum(grp) over(order by first_number) as grp
    821        from
    822        (
    823           select
    824               bin AS first_number,
    825               LEAD(bin) AS last_number
    826               IF(first_number <> lag(bin, 2) over (order by bin) + 1, 1, 0) as grp
    827           from
    828               t1
    829        )
    830     )
    831     group by grp
    832     order by 1
    833     */
    834 
    835     /*
    836 WITH Hallo AS (WITH MyTest AS (WITH MyTable AS (SELECT bin FROM Test
    837 WHERE cnt>0 GROUP BY bin ORDER BY bin) SELECT bin AS first_number,
    838 LAG(bin) OVER(ORDER BY bin) AS prev, LEAD(bin) OVER(ORDER BY bin) AS
    839 last_number, IF(bin <> (LAG(bin) OVER(ORDER BY bin)+1), 1, 0) AS
    840 grp_prev, IF(bin <> (LEAD(bin) OVER(ORDER BY bin)-1), 1, 0) as grp_nxt
    841 FROM MyTable ) SELECT first_number, last_number, sum(grp_nxt)
    842 over(order by first_number) as grp FROM MyTest) SELECT * FROM Hallo
    843 
    844 
    845     WITH MyTable AS
    846 (SELECT bin FROM Test WHERE cnt>0 GROUP BY bin ORDER BY bin)
    847 SELECT bin, LAG(bin) OVER(ORDER BY bin) AS prev, LEAD(bin) OVER(ORDER BY bin) AS nxt,
    848 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
    849 
    850 "-1","0"
    851 "0","0"
    852 "2","3"       2-4
    853 "3","5"
    854 "4","6"
    855 "5","0"
    856 "6","0"
    857 "7","7"      7-7
    858 "8","0"
    859 "9","9"      9-12
    860 "10","10"
    861 "11","11"
    862 "12","12"
    863 "13","0"
    864 "14","0"
    865 "15","0"
    866 "16","16"   16-19
    867 "17","17"
    868 "18","18"
    869 "19","19"
    870     */
    871 
    872861    // -------------------------------------------------------------------
    873862    // -------------------------- MonteCarloFiles ------------------------
     
    878867    Time start3;
    879868
    880     /* 02:get-monte-carlo-files.sql */
    881869    mysqlpp::Query query3(&connection);
    882870    query3 <<
     
    885873        "   FileId INT UNSIGNED NOT NULL,\n"
    886874        "   PRIMARY KEY (FileId) USING HASH\n"
    887         ") ENGINE=Memory\n"
     875        ") ENGINE=MEMORY\n"
    888876        "AS\n"
    889877        "(\n"
     
    897885        "      factmc.RunInfoMC\n"
    898886        "   ON\n"
    899         //"      ThetaMin BETWEEN lo AND hi OR ThetaMax BETWEEN lo AND hi\n" // Includes BOTH edges
    900887        "      (ThetaMin>=lo AND ThetaMin<hi) OR (ThetaMax>lo AND ThetaMax<=hi)\n"
    901888        "   WHERE\n"
     
    905892        "      FileId\n" // In order: faster
    906893        ")";
    907 /*
    908         "   SELECT\n"
    909         "      FileId\n"
    910         "   FROM\n"
    911         "      factmc.RunInfoMC\n"
    912         "   WHERE\n"
    913         "      PartId=1 AND\n"
    914         "      FileId%%2=0 AND\n"
    915         "      (%0:range)\n"
    916         ")";
    917 */
    918894
    919895    query3.parse();
     
    941917    cout << separator("MonteCarloArea") << '\n';
    942918
    943     Time start3b;
    944 
    945     /* 02:get-monte-carlo-files.sql */
    946     mysqlpp::Query query3b(&connection);
    947     query3b <<
    948         "CREATE TEMPORARY TABLE MonteCarloArea ENGINE=Memory\n"
     919    Time start4;
     920
     921    mysqlpp::Query query4(&connection);
     922    query4 <<
     923        "CREATE TEMPORARY TABLE MonteCarloArea ENGINE=MEMORY\n"
    949924        "AS\n"
    950925        "(\n"
     
    964939        ")";
    965940
    966     query3b.parse();
     941    query4.parse();
    967942    for (auto it=env.cbegin(); it!=env.cend(); it++)
    968         query3b.template_defaults[it->first.c_str()] = it->second.c_str();
     943        query4.template_defaults[it->first.c_str()] = it->second.c_str();
    969944
    970945    if (print_queries)
    971         PrintQuery(query3b.str());
    972 
    973     qlog << query3b << ";\n" << endl;
     946        PrintQuery(query4.str());
     947
     948    qlog << query4 << ";\n" << endl;
    974949    if (connection.connected())
    975950    {
    976         cout << query3b.execute().info() << endl;
     951        cout << query4.execute().info() << endl;
    977952        ShowWarnings(connection);
    978953        if (Dump(flog, connection, "MonteCarloArea")!=1)
     
    982957        }
    983958
    984         const auto sec3b = Time().UnixTime()-start3b.UnixTime();
    985         cout << "Execution time: " << sec3b << "s\n\n";
    986     }
    987 
    988 
    989     // -------------------------------------------------------------------
    990     // ------------------------ ThetaDistribution ------------------------
    991     // -------------------------------------------------------------------
    992 
    993     cout << separator("ThetaDistribution") << '\n';
    994 
    995     Time start4;
    996 
    997     /* 02:get-theta-distribution.sql */
    998     mysqlpp::Query query4(&connection);
    999     query4 <<
    1000         "CREATE TEMPORARY TABLE ThetaDistribution\n"
     959        const auto sec4 = Time().UnixTime()-start4.UnixTime();
     960        cout << "Execution time: " << sec4 << "s\n\n";
     961    }
     962
     963
     964    // -------------------------------------------------------------------
     965    // ----------------------------- SummaryMC ---------------------------
     966    // -------------------------------------------------------------------
     967
     968    cout << separator("SummaryOriginalMC") << '\n';
     969
     970    Time start5;
     971
     972    // This table combines the analysis results vs. Binning in Estimated Energy and Simulated Energy
     973    mysqlpp::Query query5(&connection);
     974    query5 <<
     975        "CREATE TEMPORARY TABLE Summary%100:table\n"
    1001976        "(\n"
    1002         "   `.theta`    SMALLINT UNSIGNED NOT NULL,\n"
    1003         "   lo          DOUBLE            NOT NULL COMMENT 'Lower edge of zenith distance bin in degree',\n"
    1004         "   hi          DOUBLE            NOT NULL COMMENT 'Upper edge of zenith distance bin in degree',\n"
    1005         "   CountN      INT UNSIGNED      NOT NULL,\n"
    1006         "   ErrCountN   DOUBLE            NOT NULL,\n"
    1007         "   OnTime      FLOAT             NOT NULL,\n"
    1008         "   ZdWeight    DOUBLE            NOT NULL COMMENT 'tau(delta theta)',\n"
    1009         "   ErrZdWeight DOUBLE            NOT NULL COMMENT 'sigma(tau)',\n"
    1010         "   PRIMARY KEY (`.theta`) USING HASH\n"
    1011         ") ENGINE=Memory\n"
     977        "   `.theta`      SMALLINT UNSIGNED NOT NULL,\n"
     978        "   `.sparse_sim` SMALLINT UNSIGNED NOT NULL,\n"
     979        "   `.dense_sim`  SMALLINT UNSIGNED NOT NULL,\n"
     980        "   CountN        DOUBLE            NOT NULL,\n"
     981        "   SumW          DOUBLE            NOT NULL,\n"
     982        "   SumW2         DOUBLE            NOT NULL,\n"
     983        "   INDEX (`.theta`)      USING HASH,\n"
     984        "   INDEX (`.sparse_sim`) USING HASH,\n"
     985        "   INDEX (`.dense_sim`) USING HASH\n"
     986        ") ENGINE=MEMORY\n"
    1012987        "AS\n"
    1013988        "(\n"
    1014         "   WITH EventCount AS\n"
     989        "   WITH Table0 AS\n"
    1015990        "   (\n"
    1016991        "      SELECT\n"
    1017         "         INTERVAL(DEGREES(Theta), %100:bins) AS `.theta`,\n"
    1018         "         COUNT(*) AS CountN\n"
     992        "         INTERVAL(%101:column, %102:theta) AS `.theta`,\n"
     993        "         INTERVAL(LOG10(Energy), %103:sparse) AS `.sparse_sim`,\n"
     994        "         INTERVAL(LOG10(Energy), %104:dense)  AS `.dense_sim`,\n"
     995        "         (%105:spectrum)/pow(Energy, SpectralIndex) AS SpectralWeight\n"
    1019996        "      FROM\n"
    1020997        "         MonteCarloFiles\n"
    1021998        "      LEFT JOIN\n"
    1022         "         factmc.OriginalMC USING(FileId)\n"
     999        "         factmc.RunInfoMC USING (FileId)\n"
     1000        "      LEFT JOIN\n"
     1001        "         factmc.%100:table USING (FileId)\n"
     1002        "   )\n"
     1003        "   SELECT\n"
     1004        "      `.theta`,\n"
     1005        "      `.sparse_sim`,\n"
     1006        "      `.dense_sim`,\n"
     1007        "      COUNT(*)                    AS  CountN,\n"
     1008        "      SUM(    SpectralWeight   )  AS  SumW,\n"
     1009        "      SUM(POW(SpectralWeight,2))  AS  SumW2\n"
     1010        "   FROM\n"
     1011        "      Table0\n"
     1012        "   GROUP BY\n"
     1013        "      `.theta`, `.sparse_sim`, `.dense_sim`\n"
     1014        ")";
     1015
     1016    query5.parse();
     1017    for (auto it=env.cbegin(); it!=env.cend(); it++)
     1018        query5.template_defaults[it->first.c_str()] = it->second.c_str();
     1019
     1020    query5.template_defaults["table"]    = "OriginalMC";
     1021    query5.template_defaults["column"]   = "DEGREES(Theta)";
     1022    query5.template_defaults["sparse"]   = str_sparse.c_str();
     1023    query5.template_defaults["dense"]    = str_dense.c_str();
     1024    query5.template_defaults["theta"]    = str_theta.c_str();
     1025    query5.template_defaults["spectrum"] = spectrum.c_str();
     1026
     1027    if (print_queries)
     1028        PrintQuery(query5.str());
     1029
     1030    qlog << query5 << ";\n" << endl;
     1031    if (connection.connected())
     1032    {
     1033        cout << query5.execute().info() << endl;
     1034        ShowWarnings(connection);
     1035        Dump(flog, connection, "SummaryOriginalMC");
     1036
     1037        const auto sec5 = Time().UnixTime()-start5.UnixTime();
     1038        cout << "Execution time: " << sec5 << "s\n\n";
     1039    }
     1040
     1041    // -------------------------------------------------------------------
     1042
     1043    cout << separator("SummaryEventsMC") << '\n';
     1044
     1045    Time start5b;
     1046
     1047    query5.template_defaults["table"]  = "EventsMC";
     1048    query5.template_defaults["column"] = "Zd";
     1049
     1050    if (print_queries)
     1051        PrintQuery(query5.str());
     1052
     1053    qlog << query5 << ";\n" << endl;
     1054    if (connection.connected())
     1055    {
     1056        cout << query5.execute().info() << endl;
     1057        ShowWarnings(connection);
     1058        Dump(flog, connection, "SummaryEventsMC");
     1059
     1060        const auto sec5b = Time().UnixTime()-start5b.UnixTime();
     1061        cout << "Execution time: " << sec5b << "s\n\n";
     1062    }
     1063
     1064    // -------------------------------------------------------------------
     1065    // ---------------------------- ThetDist -----------------------------
     1066    // -------------------------------------------------------------------
     1067
     1068    Time start6;
     1069
     1070    // This table combines the analysis results vs. Binning in Estimated Energy and Simulated Energy
     1071    mysqlpp::Query query6(&connection);
     1072    query6 <<
     1073        "CREATE TEMPORARY TABLE ThetaDist\n"
     1074        "ENGINE=MEMORY\n"
     1075        "AS\n"
     1076        "(\n"
     1077        "   WITH ThetaCount AS\n"
     1078        "   (\n"
     1079        "      SELECT\n"
     1080        "         `.theta`,\n"
     1081        "         SUM(CountN) AS CountN\n"
     1082        "      FROM\n"
     1083        "         SummaryOriginalMC\n"
    10231084        "      GROUP BY\n"
    10241085        "         `.theta`\n"
    10251086        "   )\n"
    10261087        "   SELECT\n"
    1027         "      `.theta`, lo, hi,\n"
     1088        "      `.theta`,\n"
    10281089        "      CountN,\n"
    10291090        "      SQRT(CountN) AS ErrCountN,\n"
     
    10341095        "      ObservationTime\n"
    10351096        "   LEFT JOIN\n"
    1036         "      EventCount USING(`.theta`)\n"
     1097        "      ThetaCount USING (`.theta`)\n"
    10371098        "   LEFT JOIN\n"
    10381099        "      BinningTheta ON `.theta`=bin\n"
     
    10411102        ")";
    10421103
    1043     query4.parse();
     1104    query6.parse();
    10441105    for (auto it=env.cbegin(); it!=env.cend(); it++)
    1045         query4.template_defaults[it->first.c_str()] = it->second.c_str();
    1046 
    1047     query4.template_defaults["bins"] = strzd.c_str();
     1106        query6.template_defaults[it->first.c_str()] = it->second.c_str();
    10481107
    10491108    if (print_queries)
    1050         PrintQuery(query4.str());
    1051 
    1052     qlog << query4 << ";\n" << endl;
     1109        PrintQuery(query6.str());
     1110
     1111    qlog << query6 << ";\n" << endl;
    10531112    if (connection.connected())
    10541113    {
    1055         cout << query4.execute().info() << endl;
     1114        cout << query6.execute().info() << endl;
    10561115        ShowWarnings(connection);
    1057         Dump(flog, connection, "ThetaDistribution");
    1058 
    1059         const auto sec4 = Time().UnixTime()-start4.UnixTime();
    1060         cout << "Execution time: " << sec4 << "s\n";
    1061 
    1062 
    1063         if (verbose>0)
    1064         {
    1065             const mysqlpp::StoreQueryResult res4 = connection.query("SELECT * FROM ThetaDistribution").store();
    1066 
    1067             cout << "  Bin   MC Counts     OnTime ZdWeight\n";
    1068             const auto bins = binning_theta.vec();
    1069             for (auto ir=res4.cbegin(); ir!=res4.cend(); ir++)
    1070             {
    1071                 const mysqlpp::Row &row = *ir;
    1072 
    1073                 const uint32_t bin = row[".theta"];
    1074                 cout << setw(5) << bins[bin] << ": " << setw(10) << row["CountN"] << " " << setw(10) << row["OnTime"] << " " << setw(10) << row["ZdWeight"] << '\n';
    1075             }
    1076             cout << endl;
    1077         }
     1116        Dump(flog, connection, "ThetaDist");
     1117
     1118        const auto sec6 = Time().UnixTime()-start6.UnixTime();
     1119        cout << "Execution time: " << sec6 << "s\n\n";
    10781120    }
    10791121
     
    10841126             .binningx = binning_theta,
    10851127             .binningy = {},
    1086              .table    = "ThetaDistribution",
     1128             .table    = "ThetaDist",
    10871129             .x        = ".theta",
    10881130             .y        = "",
     
    11011143             .binningx = binning_theta,
    11021144             .binningy = {},
    1103              .table    = "ThetaDistribution",
     1145             .table    = "ThetaDist",
    11041146             .x        = ".theta",
    11051147             .y        = "",
     
    11181160             .binningx = binning_theta,
    11191161             .binningy = {},
    1120              .table    = "ThetaDistribution",
     1162             .table    = "ThetaDist",
    11211163             .x        = ".theta",
    11221164             .y        = "",
     
    11291171         });
    11301172
    1131 
    1132     // -------------------------------------------------------------------
    1133     // ------------------------- 5: AnalysisData -------------------------
    1134     // -------------------------------------------------------------------
    1135 
    1136     cout << separator("AnalysisData") << '\n';
    1137 
    1138     Time start5;
    1139 
    1140     /* 02:analyze-data.sql */
    1141     mysqlpp::Query query5(&connection);
    1142     sindent indent5(query5);
    1143     query5 <<
    1144         "CREATE TEMPORARY TABLE AnalysisData\n"
     1173    // -------------------------------------------------------------------
     1174    // --------------------------- WeightedMC ----------------------------
     1175    // -------------------------------------------------------------------
     1176
     1177    Time start7;
     1178
     1179    // This table combines the analysis results vs. Binning in Estimated Energy and Simulated Energy
     1180    mysqlpp::Query query7(&connection);
     1181    query7 <<
     1182        "CREATE TEMPORARY TABLE Weighted%100:table (\n"
     1183        "      `.theta`      SMALLINT UNSIGNED NOT NULL,\n"
     1184        "      `.sparse_sim` SMALLINT UNSIGNED NOT NULL,\n"
     1185        "      `.dense_sim`  SMALLINT UNSIGNED NOT NULL,\n"
     1186        "      CountN        DOUBLE            NOT NULL,\n"
     1187        "      SumW          DOUBLE            NOT NULL,\n"
     1188        "      SumW2         DOUBLE            NOT NULL,\n"
     1189        "      INDEX (`.theta`)      USING HASH,\n"
     1190        "      INDEX (`.sparse_sim`) USING HASH,\n"
     1191        "      INDEX (`.dense_sim`)  USING HASH\n"
     1192        ")\n"
     1193        "ENGINE=MEMORY\n"
     1194        "AS\n"
    11451195        "(\n"
    1146         "   `.energy`       SMALLINT UNSIGNED NOT NULL,\n"
    1147         "   `Signal`        DOUBLE            NOT NULL,\n"
    1148         "   `Background`    DOUBLE            NOT NULL,\n"
    1149         "   `Excess`        DOUBLE            NOT NULL,\n"
    1150         "   `Significance`  DOUBLE            NOT NULL,\n"
    1151         "   `ErrExcess`     DOUBLE            NOT NULL,\n"
     1196        "   SELECT\n"
     1197        "      `.theta`,\n"
     1198        "      `.sparse_sim`,\n"
     1199        "      `.dense_sim`,\n"
     1200        "      S.CountN,\n"
     1201        "      S.SumW*ZdWeight AS SumW,\n"
     1202        "      S.SumW2*POW(ErrZdWeight, 2) AS SumW2\n"
     1203        "   FROM\n"
     1204        "      Summary%100:table S\n"
     1205        "   INNER JOIN\n"
     1206        "      ThetaDist USING (`.theta`)\n"
     1207        ")";
     1208
     1209    query7.parse();
     1210    for (auto it=env.cbegin(); it!=env.cend(); it++)
     1211        query7.template_defaults[it->first.c_str()] = it->second.c_str();
     1212
     1213    query7.template_defaults["table"] = "OriginalMC";
     1214
     1215    if (print_queries)
     1216        PrintQuery(query7.str());
     1217
     1218    qlog << query7 << ";\n" << endl;
     1219    if (connection.connected())
     1220    {
     1221        cout << query7.execute().info() << endl;
     1222        ShowWarnings(connection);
     1223        Dump(flog, connection, "WeightedOriginalMC");
     1224
     1225        const auto sec7 = Time().UnixTime()-start7.UnixTime();
     1226        cout << "Execution time: " << sec7 << "s\n\n";
     1227    }
     1228
     1229    // -------------------------------------------------------------------
     1230
     1231    Time start7b;
     1232
     1233    query7.template_defaults["table"] = "EventsMC";
     1234
     1235    if (print_queries)
     1236        PrintQuery(query7.str());
     1237
     1238    qlog << query7 << ";\n" << endl;
     1239    if (connection.connected())
     1240    {
     1241        cout << query7.execute().info() << endl;
     1242        ShowWarnings(connection);
     1243        Dump(flog, connection, "WeightedEventsMC");
     1244
     1245        const auto sec7b = Time().UnixTime()-start7b.UnixTime();
     1246        cout << "Execution time: " << sec7b << "s\n\n";
     1247    }
     1248
     1249    // -------------------------------------------------------------------
     1250    // --------------------------- AnalysisMC ----------------------------
     1251    // -------------------------------------------------------------------
     1252
     1253    cout << separator("AnalysisMC") << '\n';
     1254
     1255    Time start8;
     1256
     1257    // This table combines the analysis results vs. Binning in Estimated Energy and Simulated Energy
     1258    mysqlpp::Query query8(&connection);
     1259    sindent indent8(query8);
     1260    query8 <<
     1261/*
     1262        "CREATE TEMPORARY TABLE AnalysisMC\n"
     1263        "(\n"
     1264        "   `.sparse`     SMALLINT UNSIGNED NOT NULL,\n"
     1265        "   `.dense`      SMALLINT UNSIGNED NOT NULL,\n"
     1266        "   SignalW       DOUBLE            NOT NULL,\n" // vs Eest (for spectral analysis)
     1267        "   SignalN       DOUBLE            NOT NULL,\n" // vs Eest
     1268        "   BackgroundN   DOUBLE            NOT NULL,\n" // vs Eest
     1269        "   BackgroundW   DOUBLE            NOT NULL,\n" // vs Eest
     1270        "   ExcessW       DOUBLE            NOT NULL,\n" // vs Eest
     1271        "   ExcessN       DOUBLE            NOT NULL,\n" // vs Eest
     1272        "   ErrExcessW    DOUBLE            NOT NULL,\n" // vs Eest
     1273        "   ErrExcessN    DOUBLE            NOT NULL,\n" // vs Eest
     1274        "   ThresholdW    DOUBLE            NOT NULL,\n" // vs Esim (for simulation analysis)
     1275        "   ThresholdN    DOUBLE            NOT NULL,\n" // vs Esim
     1276        "   BiasEst       DOUBLE            NOT NULL,\n" // vs Eest
     1277        "   BiasSim       DOUBLE            NOT NULL,\n" // vs Esim
     1278        "   ResolutionEst DOUBLE,\n"
     1279        "   ResolutionSim DOUBLE,\n"
     1280        "   INDEX (`.sparse`) USING HASH,\n"
     1281        "   INDEX (`.dense`)  USING HASH\n"
     1282        ") ENGINE=Memory\n"
     1283*/
     1284        "CREATE TEMPORARY TABLE AnalysisMC ENGINE=MEMORY\n"
     1285        "AS\n"
     1286        "(\n"
     1287        "   WITH Excess AS\n"
     1288        "   (\n"                            << indent(6)
     1289        << ifstream(analysis_sql).rdbuf()   << indent(0) <<
     1290        "   ),\n"                           <<
     1291        "   Result AS\n"
     1292        "   (\n"                            << indent(6)
     1293        << simulation_sql << indent(0) << // Must end with EOL and not in the middle of a comment
     1294        "   )\n"
     1295        "   SELECT * FROM Result\n"
     1296        ")";
     1297
     1298    query8.parse();
     1299    for (auto it=env.cbegin(); it!=env.cend(); it++)
     1300        query8.template_defaults[it->first.c_str()] = it->second.c_str();
     1301
     1302    //query6.template_defaults["columns"]   = "FileId, EvtNumber, CorsikaNumReuse,";
     1303    query8.template_defaults["columns"]   = "Zd, Energy, SpectralIndex,";
     1304    query8.template_defaults["files"]     = "MonteCarloFiles";
     1305    query8.template_defaults["runinfo"]   = "factmc.RunInfoMC";
     1306    query8.template_defaults["events"]    = "factmc.EventsMC";
     1307    query8.template_defaults["positions"] = "factmc.PositionMC";
     1308
     1309    query8.template_defaults["sparse"]    = str_sparse.c_str();
     1310    query8.template_defaults["dense"]     = str_dense.c_str();
     1311    query8.template_defaults["theta"]     = str_theta.c_str();
     1312    query8.template_defaults["spectrum"]  = spectrum.c_str();
     1313    query8.template_defaults["estimator"] = estimator.c_str();
     1314
     1315    if (print_queries)
     1316        PrintQuery(query8.str());
     1317
     1318    qlog << query8 << ";\n" << endl;
     1319    if (connection.connected())
     1320    {
     1321        cout << query8.execute().info() << endl;
     1322        ShowWarnings(connection);
     1323        Dump(flog, connection, "AnalysisMC");
     1324
     1325        const auto sec8 = Time().UnixTime()-start8.UnixTime();
     1326        cout << "Execution time: " << sec8 << "s\n\n";
     1327    }
     1328
     1329
     1330    // -------------------------------------------------------------------
     1331    // ------------------------- SimulatedSpectrum -----------------------
     1332    // -------------------------------------------------------------------
     1333
     1334    const vector<string> binnings = { "dense", "sparse", "theta" };
     1335
     1336    for (auto ib=binnings.cbegin(); ib!=binnings.cend(); ib++)
     1337    {
     1338        const string table = "Summary"+(*ib=="theta" ? "Theta" : "TrueEnergy_"+*ib);
     1339
     1340        cout << separator(table) << '\n';
     1341
     1342        // [lo^(1+gammaa) - hi^(1+gamma)] / (1+gamma)
     1343
     1344        Time start9;
     1345
     1346        /*
     1347        "CREATE TEMPORARY TABLE SummarySimulated\n"
     1348        "(\n"
     1349        "   `.energy` SMALLINT UNSIGNED NOT NULL COMMENT 'Bin Index [MC Energy]',\n"
     1350        "   CountN    DOUBLE            NOT NULL,\n"  // COMMENT 'Number of events',\n"
     1351        "   CountW    DOUBLE            NOT NULL,\n"  // COMMENT 'Sum of weights, reweighted sum',\n"
     1352        "   CountW2   DOUBLE            NOT NULL,\n"  // COMMENT 'Sum of square of weights'\n"
    11521353        "   PRIMARY KEY (`.energy`) USING HASH\n"
    11531354        ") ENGINE=Memory\n"
     1355        */
     1356
     1357        mysqlpp::Query query9(&connection);
     1358        sindent indent9(query9);
     1359        query9 <<
     1360            "CREATE TEMPORARY TABLE %100:table ENGINE=MEMORY "
     1361            "AS\n"
     1362            "(\n"
     1363            << indent(3) << summary_sim_sql << indent(0) <<
     1364            ")";
     1365
     1366        // [ Sa^2/a^2 + Sb^2/b^2 ] * a^2/b^2
     1367        // [ (sc^2)/c^2+(sb^2)/b^2+(sa^2)/a^2 ]  * a^2*b^2/c^2
     1368
     1369        // [ a/b - c^2/d^2] --> (4*Sc^2*c^2)/d^4 + (4*Sd^2*c^4)/d^6 + Sa^2/b^2 + (Sb^2*a^2)/b^4
     1370
     1371        query9.parse();
     1372        for (auto it=env.cbegin(); it!=env.cend(); it++)
     1373            query9.template_defaults[it->first.c_str()] = it->second.c_str();
     1374
     1375        query9.template_defaults["table"]    = table.c_str();
     1376        query9.template_defaults["binning"]  = *ib=="theta" ? "BinningTheta" : ("BinningEnergy_"+*ib).c_str();
     1377        query9.template_defaults["bin"]      = *ib=="theta" ? "`.theta`"     : ("`."+*ib+"_sim`").c_str();
     1378        query9.template_defaults["binwidth"] = *ib=="theta" ? "1"            : "1000/(POW(10,hi)-POW(10,lo))";
     1379
     1380        if (print_queries)
     1381            PrintQuery(query9.str());
     1382
     1383        qlog << query9 << ";\n" << endl;
     1384        if (connection.connected())
     1385        {
     1386            cout << query9.execute().info() << endl;
     1387            ShowWarnings(connection);
     1388            Dump(flog, connection, table);
     1389
     1390            const auto sec9 = Time().UnixTime()-start9.UnixTime();
     1391            cout << "Execution time: " << sec9 << "s\n";
     1392        }
     1393
     1394        Histogram hist_sim;
     1395        hist_sim.table = table;
     1396        hist_sim.dir   = *ib=="theta" ? "MC/theta" : "MC/"+*ib+"/TrueEnergy";
     1397        hist_sim.x     = *ib=="theta" ? ".theta" : "."+*ib+"_sim";
     1398        hist_sim.axisx = *ib=="theta" ? "Zenith Angle #theta [#circ]" : "Energy lg(E_{sim}/GeV)";
     1399        hist_sim.stats = false;
     1400
     1401        if (*ib=="theta")
     1402            hist_sim.binningx=binning_theta;
     1403        if (*ib=="dense")
     1404            hist_sim.binningx=binning_dense;
     1405        if (*ib=="sparse")
     1406            hist_sim.binningx=binning_sparse;
     1407
     1408        hist_sim.axisy = "Counts";
     1409
     1410        hist_sim.title = "";
     1411        hist_sim.v     = "SimCountN";
     1412        hist_sim.err   = "ErrSimCountN";
     1413        WriteHistogram(connection, hist_sim);
     1414
     1415        hist_sim.title = "";
     1416        hist_sim.v     = "TrigCountN";
     1417        hist_sim.err   = "ErrTrigCountN";
     1418        WriteHistogram(connection, hist_sim);
     1419
     1420        hist_sim.title = "";
     1421        hist_sim.v     = "SignalN";
     1422        hist_sim.err   = "ErrSignalN";
     1423        WriteHistogram(connection, hist_sim);
     1424
     1425
     1426        hist_sim.axisy = *ib=="theta"?"dN/dE [cm^{-2} s^{-1}]":"dN/dE [cm^{-2} s^{-1} TeV^{-1}]";
     1427
     1428        hist_sim.title = "";
     1429        hist_sim.v     = "SimFluxW";
     1430        hist_sim.err   = "ErrSimFluxW";
     1431        WriteHistogram(connection, hist_sim);
     1432
     1433        hist_sim.title = "";
     1434        hist_sim.v     = "TrigFluxW";
     1435        hist_sim.err   = "ErrTrigFluxW";
     1436        WriteHistogram(connection, hist_sim);
     1437
     1438        hist_sim.title = "";
     1439        hist_sim.v     = "ExcessFluxW";
     1440        hist_sim.err   = "ErrExcessFluxW";
     1441        WriteHistogram(connection, hist_sim);
     1442
     1443        hist_sim.title = "";
     1444        hist_sim.v     = "SignalFluxW";
     1445        hist_sim.err   = "ErrSignalFluxW";
     1446        WriteHistogram(connection, hist_sim);
     1447
     1448        hist_sim.title = "";
     1449        hist_sim.v     = "BackgroundFluxW";
     1450        hist_sim.err   = "ErrBackgroundFluxW";
     1451        WriteHistogram(connection, hist_sim);
     1452
     1453
     1454        hist_sim.title = "";
     1455        hist_sim.v     = "AvgEnergySimW";
     1456        hist_sim.err   = "";
     1457        hist_sim.axisy = "<E_{sim}>/GeV";
     1458        WriteHistogram(connection, hist_sim);
     1459
     1460        hist_sim.title = "";
     1461        hist_sim.v     = "AvgEnergyEstW";
     1462        hist_sim.err   = "";
     1463        hist_sim.axisy = "<E_{sim}>/GeV";
     1464        WriteHistogram(connection, hist_sim);
     1465
     1466
     1467        hist_sim.title = "";
     1468        hist_sim.v     = "EffectiveAreaN";
     1469        hist_sim.err   = "ErrEffectiveAreaN";
     1470        hist_sim.axisy = "A_{eff} [cm^{2}]";
     1471        WriteHistogram(connection, hist_sim);
     1472
     1473        hist_sim.title = "";
     1474        hist_sim.v     = "EffectiveAreaW";
     1475        hist_sim.err   = "ErrEffectiveAreaW";
     1476        hist_sim.axisy = "A_{eff} [cm^{2}]";
     1477        WriteHistogram(connection, hist_sim);
     1478
     1479
     1480        hist_sim.title = "";
     1481        hist_sim.v     = "BiasW";
     1482        hist_sim.err   = "ErrBiasW";
     1483        hist_sim.axisy = "<lg(E_{est}/E_{sim})>";
     1484        WriteHistogram(connection, hist_sim);
     1485
     1486        hist_sim.title = "";
     1487        hist_sim.v     = "ResolutionW";
     1488        hist_sim.err   = "";
     1489        hist_sim.axisy = "#sigma(lg(E_{est}/E_{sim}))";
     1490        WriteHistogram(connection, hist_sim);
     1491
     1492        hist_sim.title = "";
     1493        hist_sim.v     = "TriggerEfficiencyN";
     1494        hist_sim.err   = "ErrTriggerEfficiencyN";
     1495        hist_sim.axisy = "Efficiency";
     1496        WriteHistogram(connection, hist_sim);
     1497
     1498        hist_sim.title = "";
     1499        hist_sim.v     = "TriggerEfficiencyW";
     1500        hist_sim.err   = "ErrTriggerEfficiencyW";
     1501        hist_sim.axisy = "Efficiency";
     1502        WriteHistogram(connection, hist_sim);
     1503
     1504        hist_sim.title = "";
     1505        hist_sim.v     = "CutEfficiencyN";
     1506        hist_sim.err   = "ErrCutEfficiencyN";
     1507        hist_sim.axisy = "Efficiency";
     1508        WriteHistogram(connection, hist_sim);
     1509
     1510        hist_sim.title = "";
     1511        hist_sim.v     = "CutEfficiencyW";
     1512        hist_sim.err   = "ErrCutEfficiencyW";
     1513        hist_sim.axisy = "Efficiency";
     1514        WriteHistogram(connection, hist_sim);
     1515
     1516        if (*ib=="theta")
     1517            continue;
     1518
     1519        // -------------------------------------------------------------------
     1520        // ------------------------- SimulatedSpectrum -----------------------
     1521        // -------------------------------------------------------------------
     1522
     1523        cout << separator("SummaryEstimatedEnergy_"+*ib) << '\n';
     1524
     1525        // [lo^(1+gammaa) - hi^(1+gamma)] / (1+gamma)
     1526
     1527        Time start10;
     1528
     1529        mysqlpp::Query query10(&connection);
     1530        sindent indent10(query10);
     1531        query10 <<
     1532            "CREATE TEMPORARY TABLE SummaryEstimatedEnergy_%100:binning ENGINE=MEMORY "
     1533            "AS\n"
     1534            "(\n"
     1535            << indent(3) << summary_est_sql << indent(6) <<
     1536            ")";
     1537
     1538        query10.parse();
     1539        for (auto it=env.cbegin(); it!=env.cend(); it++)
     1540            query10.template_defaults[it->first.c_str()] = it->second.c_str();
     1541
     1542        query10.template_defaults["binning"] = ib->c_str();
     1543
     1544        if (print_queries)
     1545            PrintQuery(query10.str());
     1546
     1547        qlog << query10 << ";\n" << endl;
     1548        if (connection.connected())
     1549        {
     1550            cout << query10.execute().info() << endl;
     1551            ShowWarnings(connection);
     1552            Dump(flog, connection, "SummaryEstimatedEnergy_"+*ib);
     1553
     1554            const auto sec10 = Time().UnixTime()-start10.UnixTime();
     1555            cout << "Execution time: " << sec10 << "s\n";
     1556        }
     1557
     1558        Histogram hist_est;
     1559        hist_est.dir      = "MC/"+*ib+"/EstimatedEnergy";
     1560        hist_est.binningx = *ib=="dense"?binning_dense:binning_sparse;
     1561        hist_est.table    = "SummaryEstimatedEnergy_"+*ib;
     1562        hist_est.x        = "."+*ib+"_est";
     1563        hist_est.axisx    = "Energy lg(E_{est}/GeV)";
     1564        hist_est.stats    = false;
     1565
     1566        hist_est.axisy   = "Counts";
     1567
     1568        hist_est.title   = "";
     1569        hist_est.v       = "SignalN";
     1570        hist_est.err     = "ErrSignalN";
     1571        WriteHistogram(connection, hist_est);
     1572
     1573        hist_est.title   = "";
     1574        hist_est.v       = "BackgroundN";
     1575        hist_est.err     = "ErrBackgroundN";
     1576        WriteHistogram(connection, hist_est);
     1577
     1578        hist_est.title   = "";
     1579        hist_est.v       = "ExcessN";
     1580        hist_est.err     = "ErrExcessN";
     1581        WriteHistogram(connection, hist_est);
     1582
     1583
     1584        hist_est.title   = "";
     1585        hist_est.v       = "AvgEnergySimW";
     1586        hist_est.err     = "";
     1587        hist_est.axisy   = "<E_{sim}>/GeV";
     1588        WriteHistogram(connection, hist_est);
     1589
     1590        hist_est.title   = "";
     1591        hist_est.v       = "AvgEnergyEstW";
     1592        hist_est.err     = "";
     1593        hist_est.axisy   = "<E_{est}>/GeV";
     1594        WriteHistogram(connection, hist_est);
     1595
     1596
     1597        hist_est.axisy   = "dN/dE [cm^{-2} s^{-1} TeV^{-1}]";
     1598
     1599        hist_est.title   = "";
     1600        hist_est.v       = "SignalFluxW";
     1601        hist_est.err     = "ErrSignalFluxW";
     1602        WriteHistogram(connection, hist_est);
     1603
     1604        hist_est.title   = "";
     1605        hist_est.v       = "BackgroundFluxW";
     1606        hist_est.err     = "ErrBackgroundFluxW";
     1607        WriteHistogram(connection, hist_est);
     1608
     1609        hist_est.title   = "";
     1610        hist_est.v       = "ExcessFluxW";
     1611        hist_est.err     = "ErrExcessFluxW";
     1612        WriteHistogram(connection, hist_est);
     1613
     1614
     1615        hist_est.title   = "";
     1616        hist_est.v       = "BiasW";
     1617        hist_est.err     = "ErrBiasW";
     1618        hist_est.axisy   = "<lg(E_{est}/E_{sim})>";
     1619        WriteHistogram(connection, hist_est);
     1620
     1621        hist_est.title   = "";
     1622        hist_est.v       = "ResolutionW";
     1623        hist_est.err     = "";
     1624        hist_est.axisy   = "#sigma(lg(E_{est}/E_{sim}))";
     1625        WriteHistogram(connection, hist_est);
     1626
     1627
     1628        // -------------------------------------------------------------------
     1629        // ------------------------- SimulatedSpectrum -----------------------
     1630        // -------------------------------------------------------------------
     1631
     1632        cout << separator("EnergyMigration_"+*ib) << '\n';
     1633
     1634        Time start11;
     1635
     1636        mysqlpp::Query query11(&connection);
     1637        query11 <<
     1638            "CREATE TEMPORARY TABLE EnergyMigration_%100:binning ENGINE=MEMORY\n"
     1639            "AS\n"
     1640            "(\n"
     1641            "   SELECT\n"
     1642            "      `.%100:binning:_est`,\n"
     1643            "      `.%100:binning:_sim`,\n"
     1644            "      SUM(SignalN)  AS SignalN\n"
     1645            "   FROM\n"
     1646            "      AnalysisMC\n"
     1647            "   GROUP BY\n"
     1648            "      `.%100:binning:_est`, `.%100:binning:_sim`\n"
     1649            "   ORDER BY\n"
     1650            "      `.%100:binning:_est`, `.%100:binning:_sim`\n"
     1651        ")";
     1652
     1653        query11.parse();
     1654        for (auto it=env.cbegin(); it!=env.cend(); it++)
     1655            query11.template_defaults[it->first.c_str()] = it->second.c_str();
     1656
     1657        query11.template_defaults["binning"] = ib->c_str();
     1658
     1659        if (print_queries)
     1660            PrintQuery(query11.str());
     1661
     1662        qlog << query11 << ";\n" << endl;
     1663        if (connection.connected())
     1664        {
     1665            cout << query11.execute().info() << endl;
     1666            ShowWarnings(connection);
     1667            Dump(flog, connection, "EnergyMigration_"+*ib);
     1668
     1669            const auto sec11 = Time().UnixTime()-start11.UnixTime();
     1670            cout << "Execution time: " << sec11 << "s\n";
     1671        }
     1672
     1673        WriteHistogram(connection, {
     1674             .name     = "Migration",
     1675             .title    = "",
     1676             .dir      = "MC/"+*ib,
     1677             .binningx = *ib=="dense"?binning_dense:binning_sparse,
     1678             .binningy = *ib=="dense"?binning_dense:binning_sparse,
     1679             .table    = "EnergyMigration_"+*ib,
     1680             .x        = "."+*ib+"_sim",
     1681             .y        = "."+*ib+"_est",
     1682             .v        = "SignalN",
     1683             .err      = "",
     1684             .axisx    = "Energy lg(E_{sim}/GeV)",
     1685             .axisy    = "Energy lg(E_{est}/GeV)",
     1686             .axisz    = "Counts",
     1687             .stats    = false
     1688         });
     1689    }
     1690
     1691    if (mc_only)
     1692    {
     1693        cout << separator("Summary") << '\n';
     1694        const auto sec = Time().UnixTime()-start.UnixTime();
     1695        cout << "Total execution time: " << sec << "s\n" << endl;
     1696
     1697        return 0;
     1698    }
     1699
     1700    // -------------------------------------------------------------------
     1701    // --------------------------- SummaryData ---------------------------
     1702    // -------------------------------------------------------------------
     1703
     1704    cout << separator("SummaryData") << '\n';
     1705
     1706    Time start12;
     1707
     1708    mysqlpp::Query query12(&connection);
     1709    sindent indent12(query12);
     1710    query12 <<
     1711        "CREATE TEMPORARY TABLE SummaryData ENGINE=MEMORY\n"
     1712/*        "(\n"
     1713        "   `.theta`        SMALLINT UNSIGNED NOT NULL,\n"
     1714        "   `.sparse_est`   SMALLINT UNSIGNED NOT NULL,\n"
     1715        "   `Signal`        DOUBLE            NOT NULL,\n"
     1716        "   `ErrSignal`     DOUBLE            NOT NULL,\n"
     1717        "   `Background`    DOUBLE            NOT NULL,\n"
     1718        "   `ErrBackground` DOUBLE            NOT NULL,\n"
     1719        "   `Excess`        DOUBLE            NOT NULL,\n"
     1720        "   `ErrExcess`     DOUBLE            NOT NULL,\n"
     1721        "   `Significance`  DOUBLE            NOT NULL,\n"
     1722        "   PRIMARY KEY (`.sparse_est`) USING HASH\n"
     1723        ") ENGINE=Memory\n"*/
    11541724        "AS\n"
    11551725        "(\n"
     
    11601730        "   Result AS\n"
    11611731        "   (\n"                          << indent(6)
    1162         << ifstream(data_sql).rdbuf()     << indent(0) << // Must end with EOL and not in the middle of a comment
     1732        << data_sql << indent(0)          << // Must end with EOL and not in the middle of a comment
    11631733        "   )\n"
    11641734        "   SELECT * FROM Result\n"
    11651735        ")";
    11661736
    1167     query5.parse();
     1737    query12.parse();
    11681738    for (auto it=env.cbegin(); it!=env.cend(); it++)
    1169         query5.template_defaults[it->first.c_str()] = it->second.c_str();
     1739        query12.template_defaults[it->first.c_str()] = it->second.c_str();
    11701740
    11711741    //query5.template_defaults["columns"]   = "FileId, EvtNumber,";
    1172     query5.template_defaults["columns"]   = "";
    1173     query5.template_defaults["files"]     = "DataFiles";
    1174     query5.template_defaults["runinfo"]   = "factdata.RunInfo";
    1175     query5.template_defaults["events"]    = "factdata.Images";
    1176     query5.template_defaults["positions"] = "factdata.Position";
    1177 
    1178     query5.template_defaults["bins"]      = stre.c_str();
    1179     query5.template_defaults["estimator"] = estimator.c_str();
     1742    query12.template_defaults["columns"]   = "fZenithDistanceMean,";
     1743    query12.template_defaults["files"]     = "DataFiles";
     1744    query12.template_defaults["runinfo"]   = "factdata.RunInfo";
     1745    query12.template_defaults["events"]    = "factdata.Images";
     1746    query12.template_defaults["positions"] = "factdata.Position";
     1747
     1748    query12.template_defaults["sparse"]    = str_sparse.c_str();
     1749    query12.template_defaults["theta"]     = str_theta.c_str();
     1750    query12.template_defaults["estimator"] = estimator.c_str();
    11801751
    11811752    if (print_queries)
    1182         PrintQuery(query5.str());
    1183 
    1184     qlog << query5 << ";\n" << endl;
     1753        PrintQuery(query12.str());
     1754
     1755    qlog << query12 << ";\n" << endl;
    11851756    if (connection.connected())
    11861757    {
    1187         cout << query5.execute().info() << endl;
     1758        cout << query12.execute().info() << endl;
    11881759        ShowWarnings(connection);
    1189         Dump(flog, connection, "AnalysisData");
    1190 
    1191         const auto sec5 = Time().UnixTime()-start5.UnixTime();
    1192         cout << "Execution time: " << sec5 << "s\n";
    1193 
    1194         if (verbose>0)
    1195         {
    1196             const mysqlpp::StoreQueryResult res5 = connection.query("SELECT * FROM AnalysisData").store();
    1197 
    1198             cout << "       Bin     Signal   Background   Excess   Significance   Error" << endl;
    1199             for (auto row=res5.cbegin(); row!=res5.cend(); row++)
    1200             {
    1201                 for (auto it=row->begin(); it!=row->end(); it++)
    1202                     cout << setw(10) << *it << " ";
    1203                 cout << '\n';
    1204             }
    1205             cout << endl;
    1206         }
    1207     }
    1208 
    1209 
    1210     // -------------------------------------------------------------------
    1211     // ----------------------------- Triggered ---------------------------
    1212     // -------------------------------------------------------------------
    1213 
    1214     cout << separator("Triggered") << '\n';
    1215 
    1216     Time start6a;
    1217 
    1218     /* 02:analyze-simulation.sql */
    1219 
    1220     // This table combines the analysis results vs. Binning in Estimated Energy and Simulated Energy
    1221     mysqlpp::Query query6a(&connection);
    1222     query6a <<
    1223         "CREATE TEMPORARY TABLE Triggered\n"
    1224         "(\n"
    1225         "   `.energy` SMALLINT UNSIGNED NOT NULL,\n"
    1226         "   CountN    DOUBLE            NOT NULL,\n"
    1227         "   CountW    DOUBLE            NOT NULL,\n"
    1228         "   CountW2   DOUBLE            NOT NULL,\n"
    1229         "   INDEX (`.energy`) USING HASH\n"
    1230         ") ENGINE=Memory\n"
    1231         "AS\n"
    1232         "(\n"
    1233         "   WITH TriggeredEvents AS\n"
    1234         "   (\n"
    1235         "      SELECT\n"
    1236         "         INTERVAL(Zd, %100:theta)  AS `.theta`,\n"
    1237         "         INTERVAL(log10(Energy), %101:energy) AS `.energy`,\n"
    1238         "         (%102:spectrum)/pow(Energy, SpectralIndex) AS SpectralWeight\n"
    1239         "      FROM\n"
    1240         "         MonteCarloFiles\n"
    1241         "      LEFT JOIN\n"
    1242         "         factmc.RunInfoMC USING (FileId)\n"
    1243         "      LEFT JOIN\n"
    1244         "         factmc.EventsMC USING (FileId)\n"
    1245         "   )\n"
    1246         "   SELECT\n"
    1247         "      `.energy`,\n"
    1248         "      COUNT(*) AS CountN,\n"
    1249         "      SUM(ZdWeight*SpectralWeight) AS CountW,\n"
    1250         "      SUM(POW(ErrZdWeight*SpectralWeight,2)) AS `CountW2`\n"
    1251         "   FROM\n"
    1252         "      TriggeredEvents\n"
    1253         "   INNER JOIN\n"
    1254         "      ThetaDistribution USING(`.theta`)\n"
    1255         "   GROUP BY\n"
    1256         "      `.energy`\n"
    1257         ")";
    1258 
    1259     query6a.parse();
    1260     for (auto it=env.cbegin(); it!=env.cend(); it++)
    1261         query6a.template_defaults[it->first.c_str()] = it->second.c_str();
    1262 
    1263     query6a.template_defaults["energy"]   = strth.c_str();
    1264     query6a.template_defaults["theta"]    = strzd.c_str();
    1265     query6a.template_defaults["spectrum"] = spectrum.c_str();
    1266 
    1267     if (print_queries)
    1268         PrintQuery(query6a.str());
    1269 
    1270     qlog << query6a << ";\n" << endl;
    1271     if (connection.connected())
    1272     {
    1273         cout << query6a.execute().info() << endl;
    1274         ShowWarnings(connection);
    1275         Dump(flog, connection, "Triggered");
    1276 
    1277         const auto sec6a = Time().UnixTime()-start6a.UnixTime();
    1278         cout << "Execution time: " << sec6a << "s\n\n";
    1279     }
    1280 
    1281     // -------------------------------------------------------------------
    1282     // ----------------------------- ResultMC ----------------------------
    1283     // -------------------------------------------------------------------
    1284 
    1285     cout << separator("ResultMC") << '\n';
    1286 
    1287     Time start6;
    1288 
    1289     /* 02:analyze-simulation.sql */
    1290 
    1291     // This table combines the analysis results vs. Binning in Estimated Energy and Simulated Energy
    1292     mysqlpp::Query query6(&connection);
    1293     sindent indent6(query6);
    1294     query6 <<
    1295         "CREATE TEMPORARY TABLE AnalysisMC\n"
    1296         "(\n"
    1297         "   `.energyest`  SMALLINT UNSIGNED NOT NULL,\n"
    1298         "   `.energysim`  SMALLINT UNSIGNED NOT NULL,\n"
    1299         "   SignalW       DOUBLE            NOT NULL,\n"
    1300         "   BackgroundW   DOUBLE            NOT NULL,\n"
    1301         "   ThresholdW    DOUBLE            NOT NULL,\n"
    1302         "   SignalN       DOUBLE            NOT NULL,\n"
    1303         "   BackgroundN   DOUBLE            NOT NULL,\n"
    1304         "   ThresholdN    DOUBLE            NOT NULL,\n"
    1305         "   ExcessW       DOUBLE            NOT NULL,\n"
    1306         "   ExcessN       DOUBLE            NOT NULL,\n"
    1307         "   ErrExcessW    DOUBLE            NOT NULL,\n"
    1308         "   ErrExcessN    DOUBLE            NOT NULL,\n"
    1309         "   BiasEst       DOUBLE            NOT NULL,\n"
    1310         "   BiasSim       DOUBLE            NOT NULL,\n"
    1311         "   ResolutionEst DOUBLE,\n"
    1312         "   ResolutionSim DOUBLE,\n"
    1313         "   INDEX (`.energyest`) USING HASH,\n"
    1314         "   INDEX (`.energysim`) USING HASH\n"
    1315         ") ENGINE=Memory\n"
    1316         "AS\n"
    1317         "(\n"
    1318         "   WITH Excess AS\n"
    1319         "   (\n"                            << indent(6)
    1320         << ifstream(analysis_sql).rdbuf()   << indent(0) <<
    1321         "   ),\n"                           <<
    1322         "   Result AS\n"
    1323         "   (\n"                            << indent(6)
    1324         << ifstream(simulation_sql).rdbuf() << indent(0) << // Must end with EOL and not in the middle of a comment
    1325         "   )\n"
    1326         "   SELECT * FROM Result\n"
    1327         ")";
    1328 
    1329     query6.parse();
    1330     for (auto it=env.cbegin(); it!=env.cend(); it++)
    1331         query6.template_defaults[it->first.c_str()] = it->second.c_str();
    1332 
    1333     //query6.template_defaults["columns"]   = "FileId, EvtNumber, CorsikaNumReuse,";
    1334     query6.template_defaults["columns"]   = "Zd, Energy, SpectralIndex,";
    1335     query6.template_defaults["files"]     = "MonteCarloFiles";
    1336     query6.template_defaults["runinfo"]   = "factmc.RunInfoMC";
    1337     query6.template_defaults["events"]    = "factmc.EventsMC";
    1338     query6.template_defaults["positions"] = "factmc.PositionMC";
    1339 
    1340     query6.template_defaults["energyest"] = stre.c_str();
    1341     query6.template_defaults["energysim"] = strth.c_str();
    1342     query6.template_defaults["theta"]     = strzd.c_str();
    1343     query6.template_defaults["spectrum"]  = spectrum.c_str();
    1344     query6.template_defaults["estimator"] = estimator.c_str();
    1345 
    1346     if (print_queries)
    1347         PrintQuery(query6.str());
    1348 
    1349     qlog << query6 << ";\n" << endl;
    1350     if (connection.connected())
    1351     {
    1352         cout << query6.execute().info() << endl;
    1353         ShowWarnings(connection);
    1354         Dump(flog, connection, "AnalysisMC");
    1355 
    1356         const auto sec6 = Time().UnixTime()-start6.UnixTime();
    1357         cout << "Execution time: " << sec6 << "s\n\n";
    1358     }
    1359 
    1360     // -------------------------------------------------------------------
    1361     // ------------------------- SimulatedSpectrum -----------------------
    1362     // -------------------------------------------------------------------
    1363 
    1364     cout << separator("SimulatedSpectrum") << '\n';
    1365 
    1366     // [lo^(1+gammaa) - hi^(1+gamma)] / (1+gamma)
    1367 
    1368     Time start7;
    1369     /* 02:get-corsika-events.sql */
    1370 
    1371     // FIXME: Theta weights?
    1372     // FIXME: energysim binning
    1373     mysqlpp::Query query7(&connection);
    1374     query7 <<
    1375         "CREATE TEMPORARY TABLE SimulatedSpectrum\n"
    1376         "(\n"
    1377         "   `.energy` SMALLINT UNSIGNED NOT NULL COMMENT 'Bin Index [MC Energy]',\n"
    1378         "   CountN    DOUBLE            NOT NULL,\n"  // COMMENT 'Number of events',\n"
    1379         "   CountW    DOUBLE            NOT NULL,\n"  // COMMENT 'Sum of weights, reweighted sum',\n"
    1380         "   CountW2   DOUBLE            NOT NULL,\n"  // COMMENT 'Sum of square of weights'\n"
    1381         "   PRIMARY KEY (`.energy`) USING HASH\n"
    1382         ") ENGINE=Memory\n"
    1383         "AS\n"
    1384         "(\n"
    1385         "   SELECT\n"
    1386         "      INTERVAL(LOG10(Energy), %100:energyest) AS `.energy`,\n"
    1387         "      COUNT(*) AS CountN,\n"
    1388         "      SUM(    (%101:spectrum)/pow(Energy, SpectralIndex)   ) AS CountW,\n"   // Contents is: CountW
    1389         "      SUM(POW((%101:spectrum)/pow(Energy, SpectralIndex),2)) AS CountW2\n"   // Error    is: SQRT(CountW2)
    1390         "   FROM\n"
    1391         "      MonteCarloFiles\n"
    1392         "   LEFT JOIN\n"
    1393         "      factmc.RunInfoMC USING (FileId)\n"
    1394         "   LEFT JOIN\n"
    1395         "      factmc.OriginalMC USING (FileId)\n"
    1396         "   GROUP BY\n"
    1397         "      `.energy`\n"
    1398         "   ORDER BY\n"
    1399         "      `.energy`\n"
    1400         ")";
    1401 
    1402     query7.parse();
    1403     for (auto it=env.cbegin(); it!=env.cend(); it++)
    1404         query7.template_defaults[it->first.c_str()] = it->second.c_str();
    1405 
    1406     //query7.template_defaults["list"]      = listmc.c_str();
    1407     query7.template_defaults["energyest"] = stre.c_str();
    1408     query7.template_defaults["spectrum"]  = spectrum.c_str();
    1409 
    1410     if (print_queries)
    1411         PrintQuery(query7.str());
    1412 
    1413     qlog << query7 << ";\n" << endl;
    1414     if (connection.connected())
    1415     {
    1416         cout << query7.execute().info() << endl;
    1417         ShowWarnings(connection);
    1418         Dump(flog, connection, "SimulatedSpectrum");
    1419 
    1420         const auto sec7 = Time().UnixTime()-start7.UnixTime();
    1421         cout << "Execution time: " << sec7 << "s\n";
    1422 
    1423 
    1424         if (verbose>0)
    1425         {
    1426             const mysqlpp::StoreQueryResult res7 = connection.query("SELECT * FROM SimulatedSpectrum").store();
    1427 
    1428             cout << "       Bin CountW           CountW2" << endl;
    1429             const auto bins = binning_esim.vec();
    1430             for (auto ir=res7.cbegin(); ir!=res7.cend(); ir++)
    1431             {
    1432                 const mysqlpp::Row &row = *ir;
    1433 
    1434                 const uint32_t bin = row[".energy"];
    1435                 cout << setw(5) << bins[bin] << ": " << setw(10) << row["CountW"] << " " << setw(10) << row["CountW2"] << '\n';
    1436             }
    1437             cout << endl;
    1438         }
    1439     }
    1440 
    1441 
    1442     // -------------------------------------------------------------------
    1443     // ----------------------------- Spectrum ----------------------------
    1444     // -------------------------------------------------------------------
    1445 
    1446     cout << separator("Spectrum") << '\n';
    1447 
    1448     Time start8;
    1449 
    1450     /* 02:calculate-spectrum.sql */
    1451 
    1452     mysqlpp::Query query8(&connection);
    1453     query8 <<
    1454         // SELECT SUM(CountN) AS CountN, SUM(OnTime) AS OnTime FROM ThetaHist
    1455         "CREATE TEMPORARY TABLE Spectrum\n"
     1760        Dump(flog, connection, "SummaryData");
     1761
     1762        const auto sec12 = Time().UnixTime()-start12.UnixTime();
     1763        cout << "Execution time: " << sec12 << "s\n";
     1764    }
     1765
     1766    // -------------------------------------------------------------------
     1767    // --------------------------- Spectrum ------------------------------
     1768    // -------------------------------------------------------------------
     1769
     1770    const vector<string> targets = { "Theta", "Energy" };
     1771
     1772    for (auto ib=targets.cbegin(); ib!=targets.cend(); ib++)
     1773    {
     1774        const string table = "Spectrum"+*ib;
     1775
     1776        cout << separator(table) << '\n';
     1777
     1778        Time start13;
     1779
     1780        /*
     1781         "CREATE TEMPORARY TABLE Spectrum\n"
    14561782        "(\n"
    14571783        "   `.energy`      SMALLINT UNSIGNED NOT NULL COMMENT 'Bin Index [Energy]' PRIMARY KEY,\n"
     
    14731799        "   ErrSignalW     DOUBLE            NOT NULL COMMENT 'Error of weighted number of signal events in simulated data',\n"
    14741800        "   ErrBackgroundW DOUBLE            NOT NULL COMMENT 'Error of weighted number of background events in simulated data',\n"
    1475         "   Flux           DOUBLE            NOT NULL COMMENT 'dN/dA/dt [cm^-2 s-^1 TeV^-1]',\n"
    1476         "   ErrFlux        DOUBLE            NOT NULL COMMENT 'dN/dA/dt [cm^-2 s-^1 TeV^-1]',\n"
     1801        "   Flux           DOUBLE            NOT NULL COMMENT 'Measured Monte Carlo Flux dN/dA/dt [cm^-2 s-^1 TeV^-1]',\n"
     1802        "   ErrFlux        DOUBLE            NOT NULL COMMENT 'Error of measures Monte Carlo Flux dN/dA/dt [cm^-2 s-^1 TeV^-1]',\n"
     1803        //"   CountSim       DOUBLE            NOT NULL COMMENT 'Simulated Monte Carlo Events',\n"
     1804        //"   ErrCountSim    DOUBLE            NOT NULL COMMENT 'Error of Simulated Monte Carlo Events',\n"
     1805        //"   FluxSim        DOUBLE            NOT NULL COMMENT 'Simulated Monte Carlo Flux dN/dA/dt [cm^-2 s-^1 TeV^-1]',\n"
     1806        //"   ErrFluxSim     DOUBLE            NOT NULL COMMENT 'Error of Simulated Monte Carlo Flux dN/dA/dt [cm^-2 s-^1 TeV^-1]',\n"
    14771807        "   Bias           DOUBLE            NOT NULL COMMENT 'Energy Bias, average residual in lg(E)',\n"
    14781808        "   Resolution     DOUBLE            NOT NULL COMMENT 'Energy resolution, standard divation of residual in lg(E)',\n"
    1479         "   EfficiencyN    DOUBLE            NOT NULL COMMENT 'Simulated cut efficiency (weighted)',\n"
    1480         "   EfficiencyW    DOUBLE            NOT NULL COMMENT 'Simulated cut efficiency (unweighted)',\n"
    1481         "   ErrEfficiencyN DOUBLE            NOT NULL COMMENT 'Error of simulated cut efficiency (weighted)',\n"
    1482         "   ErrEfficiencyW DOUBLE            NOT NULL COMMENT 'Error of simulated cut efficiency (unweighted)'\n"
     1809        //"   EfficiencyN    DOUBLE            NOT NULL COMMENT 'Simulated cut efficiency (weighted)',\n"
     1810        //"   EfficiencyW    DOUBLE            NOT NULL COMMENT 'Simulated cut efficiency (unweighted)',\n"
     1811        //"   ErrEfficiencyN DOUBLE            NOT NULL COMMENT 'Error of simulated cut efficiency (weighted)',\n"
     1812        //"   ErrEfficiencyW DOUBLE            NOT NULL COMMENT 'Error of simulated cut efficiency (unweighted)'\n"
    14831813        ") ENGINE=Memory\n"
    1484         "AS\n"
    1485         "(\n"
    1486         "   WITH ThetaSums AS\n"
    1487         "   (\n"
    1488         "      SELECT\n"
    1489         "         SUM(CountN) AS CountSim,\n"
    1490         "         SUM(OnTime) AS ObsTime\n"
    1491         "      FROM\n"
    1492         "         ThetaDistribution\n"
    1493         "   ),\n"
    1494         "   Area AS\n"
    1495         "   (\n"
    1496         "      SELECT\n"
    1497         "         POW(MinImpactHi,2)*PI() AS Area\n"
    1498         "      FROM\n"
    1499         "         MonteCarloArea\n"
    1500         "   ),\n"
    1501         "   ResultMC AS\n" // Group AnalysisMC by EnergyEst Binning
    1502         "   (\n"
    1503         "      SELECT\n"
    1504         "         `.energyest`             AS `.energy`,\n"
    1505         "         ANY_VALUE(SignalW)       AS SignalW,\n"
    1506         "         ANY_VALUE(SignalW2)      AS SignalW2,\n"
    1507         "         ANY_VALUE(BackgroundW)   AS BackgroundW,\n"
    1508         "         ANY_VALUE(BackgroundW2)  AS BackgroundW2,\n"
    1509         "         ANY_VALUE(SignalN)       AS SignalN,\n"
    1510         "         ANY_VALUE(BackgroundN)   AS BackgroundN,\n"
    1511         "         ANY_VALUE(ExcessW)       AS ExcessW,\n"
    1512         "         ANY_VALUE(ExcessN)       AS ExcessN,\n"
    1513         "         ANY_VALUE(ErrExcessW)    AS ErrExcessW,\n"
    1514         "         ANY_VALUE(ErrExcessN)    AS ErrExcessN,\n"
    1515         "         ANY_VALUE(BiasEst)       AS Bias,\n"
    1516         "         ANY_VALUE(ResolutionEst) AS Resolution\n"
    1517         "      FROM\n"
    1518         "         AnalysisMC\n"
    1519         "      GROUP BY\n"
    1520         "         `.energy`\n"
    1521         "      ORDER BY\n"
    1522         "         `.energy`\n"
    1523         "   )\n"
    1524         "   SELECT\n"
    1525         "      `.energy`, lo, hi,\n"  // Scale for Theta-Weights
    1526         "      `Signal`, `Background`/5 AS `Background`, `Excess`, `ErrExcess`, `Significance`,\n"
    1527         "      SQRT(`Signal`)         AS ErrSignal,\n"
    1528         "      SQRT(`SignalW2`)       AS ErrSignalW,\n"
    1529         "      SQRT(`Background`)/5   AS ErrBackground,\n"
    1530         "      SQRT(`BackgroundW2`)/5 AS ErrBackgroundW,\n"
    1531         "      ExcessN, ExcessW, ErrExcessN, ErrExcessW, SignalW, BackgroundW,\n"
    1532         "      AnalysisData.Excess/ResultMC.ExcessW*SimulatedSpectrum.CountW * 1000/(POW(10,hi)-POW(10,lo)) /Area/ObsTime / CountSim*ObsTime AS Flux,\n"
    1533         "      AnalysisData.Excess/ResultMC.ExcessW*SimulatedSpectrum.CountW * 1000/(POW(10,hi)-POW(10,lo)) /Area/ObsTime / CountSim*ObsTime\n"
    1534         "         * SQRT(\n"
    1535         "             + POW(AnalysisData.ErrExcess / AnalysisData.Excess, 2)\n"
    1536         "             + POW(ResultMC.ErrExcessW    / ResultMC.ExcessW,    2)\n"
    1537         "             + SimulatedSpectrum.CountW2  / POW(SimulatedSpectrum.CountW,2)\n"
    1538         "           ) AS ErrFlux,\n"
    1539         "      Bias,\n"
    1540         "      Resolution,\n"
    1541         "      ResultMC.ExcessW/SimulatedSpectrum.CountW * CountSim/ObsTime AS EfficiencyW,\n"
    1542         "      ResultMC.ExcessN/SimulatedSpectrum.CountN AS EfficiencyN,\n"
    1543         "      ( POW(ResultMC.ErrExcessW/ResultMC.ExcessW, 2) + POW(SQRT(SimulatedSpectrum.CountW2)/SimulatedSpectrum.CountW, 2) )\n"
    1544         "         * POW(ResultMC.ExcessW/SimulatedSpectrum.CountW, 2) * CountSim/ObsTime AS ErrEfficiencyW,\n"
    1545         "      ( POW(ResultMC.ErrExcessN, 2) + POW(ResultMC.ExcessN, 2)/SimulatedSpectrum.CountN)/POW(SimulatedSpectrum.CountN, 2) AS ErrEfficiencyN,\n"
    1546         "      Area/10000*ResultMC.ExcessW/SimulatedSpectrum.CountW * CountSim/ObsTime AS EffAreaW,\n"
    1547         "      Area/10000*ResultMC.ExcessN/SimulatedSpectrum.CountN AS EffAreaN,\n"
    1548         "      Area/10000*( POW(ResultMC.ErrExcessW/ResultMC.ExcessW, 2) + POW(SQRT(SimulatedSpectrum.CountW2)/SimulatedSpectrum.CountW, 2) )\n"
    1549         "         * POW(ResultMC.ExcessW/SimulatedSpectrum.CountW, 2) * CountSim/ObsTime AS ErrEffAreaW,\n"
    1550         "      Area/10000*( POW(ResultMC.ErrExcessN, 2) + POW(ResultMC.ExcessN, 2)/SimulatedSpectrum.CountN)/POW(SimulatedSpectrum.CountN, 2) AS ErrEffAreaN\n"
    1551         "   FROM\n"
    1552         "      AnalysisData\n"
    1553         "   INNER JOIN\n"
    1554         "      ResultMC USING(`.energy`)\n"
    1555         "   INNER JOIN\n"
    1556         "      SimulatedSpectrum USING(`.energy`)\n"
    1557         "   INNER JOIN\n"
    1558         "      BinningEnergyEst ON `.energy`=bin\n"
    1559         "   CROSS JOIN\n"
    1560         "      ThetaSums, Area\n"
    1561         "   WHERE\n"
    1562         "      AnalysisData.Excess>0\n"
    1563         "   ORDER BY\n"
    1564         "      `.energy`\n"
    1565         ")"
    1566         ;
    1567 
    1568     // [ Sa^2/a^2 + Sb^2/b^2 ] * a^2/b^2
    1569     // [ (sc^2)/c^2+(sb^2)/b^2+(sa^2)/a^2 ]  * a^2*b^2/c^2
    1570 
    1571     query8.parse();
    1572     for (auto it=env.cbegin(); it!=env.cend(); it++)
    1573         query8.template_defaults[it->first.c_str()] = it->second.c_str();
    1574 
    1575     //query8.template_defaults["area"] = area;
    1576     //query8.template_defaults["ontime"] = resX[0]["OnTime"].data();
    1577     //query8.template_defaults["count"]  = resX[0]["CountN"].data();
    1578 
    1579     if (print_queries)
    1580         PrintQuery(query8.str());
    1581 
    1582     qlog << query8 << ";\n" << endl;
    1583     if (connection.connected())
    1584     {
    1585         cout << query8.execute().info() << endl;
    1586         ShowWarnings(connection);
    1587         Dump(flog, connection, "Spectrum");
    1588 
    1589         const auto sec8 = Time().UnixTime()-start8.UnixTime();
    1590         cout << "Execution time: " << sec8 << "s\n";
    1591 
    1592 
    1593         const mysqlpp::StoreQueryResult res8 = connection.query("SELECT * FROM Spectrum").store();
    1594 
    1595         if (verbose>0)
     1814*/
     1815
     1816        mysqlpp::Query query13(&connection);
     1817        sindent indent13(query13);
     1818        query13 <<
     1819            "CREATE TEMPORARY TABLE %100:table ENGINE=MEMORY "
     1820            "AS\n"
     1821            "(\n"
     1822            << indent(3) << spectrum_sql << indent(0) <<
     1823            ")";
     1824*/
     1825        // [ Sa^2/a^2 + Sb^2/b^2 ] * a^2/b^2
     1826        // [ (sc^2)/c^2+(sb^2)/b^2+(sa^2)/a^2 ]  * a^2*b^2/c^2
     1827
     1828        query13.parse();
     1829        for (auto it=env.cbegin(); it!=env.cend(); it++)
     1830            query13.template_defaults[it->first.c_str()] = it->second.c_str();
     1831
     1832        query13.template_defaults["table"]     = table.c_str();
     1833        query13.template_defaults["bin"]       = *ib=="Theta"       ? "`.theta`"    : "`.sparse_est`";
     1834        query13.template_defaults["id"]        = *ib=="Theta"       ? "Sim"         : "Est";
     1835        query13.template_defaults["weight"]    = *ib=="Theta"       ? "ZdWeight"    : "1";
     1836        query13.template_defaults["errweight"] = *ib=="Theta"       ? "ErrZdWeight" : "1";
     1837        if (*ib=="Theta")
    15961838        {
    1597             cout << "  Bin  Flux                   Error" << endl;
    1598             const auto bins = binning_eest.vec();
    1599             for (auto ir=res8.cbegin(); ir!=res8.cend(); ir++)
     1839            query13.template_defaults["join1"] = "SummaryTheta Sim USING (`.theta`)";
     1840            query13.template_defaults["join2"] = "ThetaDist USING (`.theta`)";
     1841        }
     1842        else
     1843        {
     1844            query13.template_defaults["join1"] = "SummaryEstimatedEnergy_sparse Est USING (`.sparse_est`)";
     1845            query13.template_defaults["join2"] = "SummaryTrueEnergy_sparse Sim ON Est.`.sparse_est`=Sim.`.sparse_sim`";
     1846        }
     1847
     1848        if (print_queries)
     1849            PrintQuery(query13.str());
     1850
     1851        qlog << query13 << ";\n" << endl;
     1852        if (connection.connected())
     1853        {
     1854            cout << query13.execute().info() << endl;
     1855            ShowWarnings(connection);
     1856            Dump(flog, connection, table);
     1857
     1858            const auto sec13 = Time().UnixTime()-start13.UnixTime();
     1859            cout << "Execution time: " << sec13 << "s\n";
     1860
     1861
     1862            const mysqlpp::StoreQueryResult res13 = connection.query("SELECT * FROM "+table).store();
     1863
     1864            // --------------------------------------------------------------------------
     1865#ifdef HAVE_ROOT
     1866            TFeldmanCousins fc;
     1867            fc.SetCL(confidence);
     1868            fc.SetMuStep(0.2);
     1869            // f.SetMuMax(10*(sig+bg)); //has to be higher than sig+bg!!
     1870            // Double_t ul=f.CalculateUpperLimit(x,y);
     1871            // x=Dat.Signal       : number of observed events in the experiment
     1872            // y=Dat.Background/5 : average number of observed events in background region
     1873
     1874            TRolke rolke(confidence);
     1875            // rolke.SetPoissonBkgBinomEff(x,y,z,tau,m)
     1876            // x=Dat.Signal     : number of observed events in the experiment
     1877            // y=Dat.Background : number of observed events in background region
     1878            // tau=0.2          : the ratio between signal and background region
     1879            // m=Sim.SimFluxN   : number of MC events generated
     1880            // z=Sim.AnaFluxN   : number of MC events observed
     1881#endif
     1882            // --------------------------------------------------------------------------
     1883
     1884            // Crab Nebula: 1 TeV: 3e-7  /  m^2 / s / TeV
     1885            // Crab Nebula: 1 TeV: 3e-11 / cm^2 / s / TeV
     1886
     1887            map<size_t, double> feldman_ul;
     1888            map<size_t, double> rolke_ul;
     1889            map<size_t, double> rolke_ll;
     1890
     1891            for (auto ir=res13.cbegin(); ir!=res13.cend(); ir++)
    16001892            {
     1893                // This is not efficient but easier to read and safer
    16011894                const mysqlpp::Row &row = *ir;
    16021895
    1603                 const uint32_t bin = row[".energy"];
    1604                 cout << setw(5) << bins[bin] << ": " << setw(10) << row["Flux"] << " " << setw(10) << row["ErrFlux"] << '\n';
     1896                const size_t bin     = row[*ib=="Theta" ? ".theta" : ".sparse_est"];
     1897                const double flux    = row["Flux"];
     1898                const double error   = row["ErrFlux"];
     1899
     1900#ifdef HAVE_ROOT
     1901                const double dat_sig = row["Signal"];
     1902                const double dat_bg  = row["Background"];
     1903
     1904                const double eff     = row["Efficiency"];
     1905                const double scale   = row["Scale"];
     1906
     1907                fc.SetMuMax(10*(dat_sig+dat_bg/5)); //has to be higher than sig+bg!!
     1908                rolke.SetPoissonBkgKnownEff(dat_sig, dat_bg, 5, eff);
     1909
     1910                if (feldman)
     1911                    feldman_ul[bin] = fc.CalculateUpperLimit(dat_sig, dat_bg/5)*scale/eff;
     1912
     1913                rolke_ll[bin] = rolke.GetLowerLimit()*scale;
     1914                rolke_ul[bin] = rolke.GetUpperLimit()*scale;
     1915#endif
     1916                if (verbose>0)
     1917                {
     1918                    cout << setw(5)  << row["center"] << ": ";
     1919                    cout << setw(10) << row["Excess"] << " ";
     1920                    cout << setw(10) << flux << " ";
     1921                    cout << setw(10) << error << " ";
     1922                    cout << setw(10) << row["Significance"] << '\n';
     1923                }
    16051924            }
    1606             cout << endl;
     1925
     1926            Histogram hist_zd;
     1927            hist_zd.dir      = "Data/"+*ib;
     1928            hist_zd.table    = table;
     1929            hist_zd.binningx = *ib=="Theta" ? binning_theta : binning_sparse;
     1930            hist_zd.x        = *ib=="Theta" ? ".theta" : ".sparse_est";
     1931            hist_zd.axisx    = *ib=="Theta" ? "Zenith Distance #theta [#circ]" : "Energy lg(E/GeV)";
     1932            hist_zd.stats    = false;
     1933
     1934            hist_zd.title    = "";
     1935            hist_zd.v        = "Signal";
     1936            hist_zd.err      = "ErrSignal";
     1937            hist_zd.axisy    = "Counts";
     1938            WriteHistogram(connection, hist_zd);
     1939
     1940            hist_zd.title    = "";
     1941            hist_zd.v        = "Background";
     1942            hist_zd.err      = "ErrBackground";
     1943            hist_zd.axisy    = "Counts";
     1944            WriteHistogram(connection, hist_zd);
     1945
     1946            hist_zd.title    = "";
     1947            hist_zd.v        = "Excess";
     1948            hist_zd.err      = "ErrExcess";
     1949            hist_zd.axisy    = "Counts";
     1950            WriteHistogram(connection, hist_zd);
     1951
     1952            hist_zd.title    = "";
     1953            hist_zd.v        = "Significance";
     1954            hist_zd.err      = "";
     1955            hist_zd.axisy    = "#sigma";
     1956            WriteHistogram(connection, hist_zd);
     1957
     1958            hist_zd.title    = "";
     1959            hist_zd.v        = "AvgEnergyEst";
     1960            hist_zd.err      = "";
     1961            hist_zd.axisy    = "<E_{est}>/GeV";
     1962            WriteHistogram(connection, hist_zd);
     1963
     1964            hist_zd.v        = "ExcessRatio";
     1965            hist_zd.err      = "ErrExcessRatio";
     1966            hist_zd.axisy    = "Ratio";
     1967            WriteHistogram(connection, hist_zd);
     1968
     1969
     1970            hist_zd.axisy    = *ib=="Theta" ? "dN/dE [cm^{-2} s^{-1}]" : "dN/dE [cm^{-2} s^{-1} TeV^{-1}]";
     1971
     1972            hist_zd.name     = "Spectrum";
     1973            hist_zd.v        = "Flux";
     1974            hist_zd.err      = "ErrFlux";
     1975            WriteHistogram(connection, hist_zd);
     1976
     1977#ifdef HAVE_ROOT
     1978            hist_zd.axisy    = *ib=="Theta" ? "UL [cm^{-2} s^{-1}]" : "UL [cm^{-2} s^{-1} TeV^{-1}]";
     1979
     1980            if (feldman)
     1981            {
     1982                hist_zd.name = "FeldmanCousins";
     1983                WriteHistogram(connection, hist_zd, feldman_ul);
     1984            }
     1985
     1986            hist_zd.name = "RolkeUL";
     1987            WriteHistogram(connection, hist_zd, rolke_ul);
     1988
     1989            hist_zd.name = "RolkeLL";
     1990            WriteHistogram(connection, hist_zd, rolke_ll);
     1991#endif
    16071992        }
    1608 
    1609         // --------------------------------------------------------------------------
    1610 
    1611         // Crab Nebula: 1 TeV: 3e-7  /  m^2 / s / TeV
    1612         // Crab Nebula: 1 TeV: 3e-11 / cm^2 / s / TeV
    1613 
    1614         sindent indentm(mlog);
    1615 
    1616         mlog << "void spectrum()\n";
    1617         mlog << "{\n" << indent(4);
    1618         mlog << "// Energy Spectrum (e.g. Crab: 3e-11 [cm^-2 s-^1 TeV^-1] @ 1TeV)\n";
    1619         mlog << "TGraphErrors g;\n";
    1620         mlog << "g.SetNameTitle(\"Spectrum\", \"Energy Spectrum\");\n";
    1621         for (auto ir=res8.cbegin(); ir!=res8.cend(); ir++)
    1622         {
    1623             // This is not efficient but easier to read and safer
    1624             const mysqlpp::Row &row = *ir;
    1625 
    1626             const double hi     = row["hi"];
    1627             const double lo     = row["lo"];
    1628             const double center = (hi+lo)/2;
    1629             const double flux   = row["Flux"];
    1630             const double error  = row["ErrFlux"];
    1631 
    1632             mlog << "g.SetPoint(g.GetN(), " << setw(8) << center << ", " << flux << ");\n";
    1633             mlog << "g.SetPointError(g.GetN()-1, 0, " << error << ");\n";
    1634         }
    1635         mlog << "g.GetXaxis()->SetTitle(\"lg(E/TeV)\");\n";
    1636         mlog << "g.GetYaxis()->SetTitle(\"dN/dE [cm^{-2} s^{-1} TeV^{-1}]\");";
    1637         mlog << endl;
    1638     }
    1639 
    1640     WriteHistogram(connection, {
    1641              .name     = "Signal",
    1642              .title    = "Signal",
    1643              .dir      = "Eest/Measurement",
    1644              .binningx = binning_eest,
    1645              .binningy = {},
    1646              .table    = "Spectrum",
    1647              .x        = ".energy",
    1648              .y        = "",
    1649              .v        = "Signal",
    1650              .err      = "ErrSignal",
    1651              .axisx    = "Energy lg(E_{est}/GeV)",
    1652              .axisy    = "Counts",
    1653              .axisz    = "",
    1654              .stats    = true
    1655          });
    1656 
    1657     WriteHistogram(connection, {
    1658              .name     = "Background",
    1659              .title    = "Background",
    1660              .dir      = "Eest/Measurement",
    1661              .binningx = binning_eest,
    1662              .binningy = {},
    1663              .table    = "Spectrum",
    1664              .x        = ".energy",
    1665              .y        = "",
    1666              .v        = "Background",
    1667              .err      = "ErrBackground",
    1668              .axisx    = "Energy lg(E_{est}/GeV)",
    1669              .axisy    = "Count",
    1670              .axisz    = "",
    1671              .stats    = true
    1672          });
    1673 
    1674     WriteHistogram(connection, {
    1675              .name     = "Excess",
    1676              .title    = "Excess",
    1677              .dir      = "Eest/Measurement",
    1678              .binningx = binning_eest,
    1679              .binningy = {},
    1680              .table    = "Spectrum",
    1681              .x        = ".energy",
    1682              .y        = "",
    1683              .v        = "Excess",
    1684              .err      = "ErrExcess",
    1685              .axisx    = "Energy lg(E_{est}/GeV)",
    1686              .axisy    = "Signal - Background (Counts)",
    1687              .axisz    = "",
    1688              .stats    = true
    1689          });
    1690 
    1691     WriteHistogram(connection, {
    1692              .name     = "SignalW",
    1693              .title    = "SignalW",
    1694              .dir      = "Eest/Simulation/Weighted",
    1695              .binningx = binning_eest,
    1696              .binningy = {},
    1697              .table    = "Spectrum",
    1698              .x        = ".energy",
    1699              .y        = "",
    1700              .v        = "SignalW",
    1701              .err      = "ErrSignalW",
    1702              .axisx    = "Energy lg(E_{est}/GeV)",
    1703              .axisy    = "Weighted",
    1704              .axisz    = "",
    1705              .stats    = true
    1706          });
    1707 
    1708     WriteHistogram(connection, {
    1709              .name     = "BackgroundW",
    1710              .title    = "BackgroundW",
    1711              .dir      = "Eest/Simulation/Weighted",
    1712              .binningx = binning_eest,
    1713              .binningy = {},
    1714              .table    = "Spectrum",
    1715              .x        = ".energy",
    1716              .y        = "",
    1717              .v        = "BackgroundW",
    1718              .err      = "ErrBackgroundW",
    1719              .axisx    = "Energy lg(E_{est}/GeV)",
    1720              .axisy    = "Weighted",
    1721              .axisz    = "",
    1722              .stats    = true
    1723          });
    1724 
    1725     WriteHistogram(connection, {
    1726              .name     = "ExcessW",
    1727              .title    = "ExcessW",
    1728              .dir      = "Eest/Simulation/Weighted",
    1729              .binningx = binning_eest,
    1730              .binningy = {},
    1731              .table    = "Spectrum",
    1732              .x        = ".energy",
    1733              .y        = "",
    1734              .v        = "ExcessW",
    1735              .err      = "ErrExcessW",
    1736              .axisx    = "Energy lg(E_{est}/GeV)",
    1737              .axisy    = "Signal - Background (Weighted)",
    1738              .axisz    = "",
    1739              .stats    = true
    1740          });
    1741 
    1742     WriteHistogram(connection, {
    1743              .name     = "Significance",
    1744              .title    = "Significance",
    1745              .dir      = "Eest/Measurement",
    1746              .binningx = binning_eest,
    1747              .binningy = {},
    1748              .table    = "Spectrum",
    1749              .x        = ".energy",
    1750              .y        = "",
    1751              .v        = "Significance",
    1752              .err      = "",
    1753              .axisx    = "Energy lg(E_{est}/GeV)",
    1754              .axisy    = "Li/Ma Significance",
    1755              .axisz    = "",
    1756              .stats    = true
    1757          });
    1758 
    1759     WriteHistogram(connection, {
    1760              .name     = "Bias",
    1761              .title    = "Energy Bias",
    1762              .dir      = "Eest",
    1763              .binningx = binning_eest,
    1764              .binningy = {},
    1765              .table    = "Spectrum",
    1766              .x        = ".energy",
    1767              .y        = "",
    1768              .v        = "Bias",
    1769              .err      = "Resolution",
    1770              .axisx    = "Energy lg(E_{sim}/GeV)",
    1771              .axisy    = "lg(E_{est}/E_{sim})",
    1772              .axisz    = "",
    1773              .stats    = false
    1774          });
    1775 
    1776     WriteHistogram(connection, {
    1777              .name     = "Resolution",
    1778              .title    = "Energy Resolution",
    1779              .dir      = "Eest",
    1780              .binningx = binning_eest,
    1781              .binningy = {},
    1782              .table    = "Spectrum",
    1783              .x        = ".energy",
    1784              .y        = "",
    1785              .v        = "Resolution",
    1786              .err      = "",
    1787              .axisx    = "Energy lg(E_{sim}/GeV)",
    1788              .axisy    = "#sigma(lg(E_{est}/E_{sim}))",
    1789              .axisz    = "",
    1790              .stats    = false
    1791          });
    1792 
    1793     WriteHistogram(connection, {
    1794              .name     = "EfficiencyN",
    1795              .title    = "Efficiency (Counts)",
    1796              .dir      = "Eest/Efficiency",
    1797              .binningx = binning_eest,
    1798              .binningy = {},
    1799              .table    = "Spectrum",
    1800              .x        = ".energy",
    1801              .y        = "",
    1802              .v        = "EfficiencyN",
    1803              .err      = "ErrEfficiencyN",
    1804              .axisx    = "Energy lg(E_{sim}/GeV)",
    1805              .axisy    = "Efficiency",
    1806              .axisz    = "",
    1807              .stats    = true
    1808          });
    1809 
    1810     WriteHistogram(connection, {
    1811              .name     = "EfficiencyW",
    1812              .title    = "Efficiency (Weighted)",
    1813              .dir      = "Eest/Efficiency",
    1814              .binningx = binning_eest,
    1815              .binningy = {},
    1816              .table    = "Spectrum",
    1817              .x        = ".energy",
    1818              .y        = "",
    1819              .v        = "EfficiencyW",
    1820              .err      = "ErrEfficiencyW",
    1821              .axisx    = "Energy lg(E_{sim}/GeV)",
    1822              .axisy    = "Efficiency",
    1823              .axisz    = "",
    1824              .stats    = true
    1825          });
    1826 
    1827     WriteHistogram(connection, {
    1828              .name     = "EffAreaN",
    1829              .title    = "Effective Area (Counts)",
    1830              .dir      = "Eest/EffArea",
    1831              .binningx = binning_eest,
    1832              .binningy = {},
    1833              .table    = "Spectrum",
    1834              .x        = ".energy",
    1835              .y        = "",
    1836              .v        = "EffAreaN",
    1837              .err      = "ErrEffAreaN",
    1838              .axisx    = "Energy lg(E_{sim}/GeV)",
    1839              .axisy    = "Effective Area A [m^{2}]",
    1840              .axisz    = "",
    1841              .stats    = false
    1842          });
    1843 
    1844     WriteHistogram(connection, {
    1845              .name     = "EffAreaW",
    1846              .title    = "Effective Area (Weighted)",
    1847              .dir      = "Eest/EffArea",
    1848              .binningx = binning_eest,
    1849              .binningy = {},
    1850              .table    = "Spectrum",
    1851              .x        = ".energy",
    1852              .y        = "",
    1853              .v        = "EffAreaW",
    1854              .err      = "ErrEffAreaW",
    1855              .axisx    = "Energy lg(E_{sim}/GeV)",
    1856              .axisy    = "Effective Area A [m^{2}]",
    1857              .axisz    = "",
    1858              .stats    = false
    1859          });
    1860 
    1861     WriteHistogram(connection, {
    1862              .name     = "Spectrum",
    1863              .title    = "Differential Energy Spectrum",
    1864              .dir      = "",
    1865              .binningx = binning_eest,
    1866              .binningy = {},
    1867              .table    = "Spectrum",
    1868              .x        = ".energy",
    1869              .y        = "",
    1870              .v        = "Flux",
    1871              .err      = "ErrFlux",
    1872              .axisx    = "Energy lg(E/GeV)",
    1873              .axisy    = "dN/dE [cm^{-2} s^{-1} TeV^{-1}]",
    1874              .axisz    = "",
    1875              .stats    = false
    1876          });
    1877 
    1878 
    1879     // -------------------------------------------------------------------
    1880     // ---------------------------- Threshold ----------------------------
    1881     // -------------------------------------------------------------------
    1882 
    1883     cout << separator("Threshold") << '\n';
    1884 
    1885     Time start9;
    1886 
    1887     /* 02:calculate-threshold.sql */
    1888 
    1889     // This query gets the analysis results versus Simulated Energy from the combined table
    1890     mysqlpp::Query query9(&connection);
    1891     query9 <<
    1892         // SELECT SUM(CountN) AS CountN, SUM(OnTime) AS OnTime FROM ThetaHist
    1893         "CREATE TEMPORARY TABLE Threshold ENGINE=Memory AS\n"
    1894         "(\n"
    1895         "   WITH ThetaSums AS\n"
    1896         "   (\n"
    1897         "      SELECT\n"
    1898         "         SUM(CountN) AS CountSim,\n"
    1899         "         SUM(OnTime) AS ObsTime\n"
    1900         "      FROM\n"
    1901         "         ThetaDistribution\n"
    1902         "   ),\n"
    1903         "   Area AS\n"
    1904         "   (\n"
    1905         "      SELECT\n"
    1906         "         POW(MinImpactHi,2)*PI() AS Area\n"
    1907         "      FROM\n"
    1908         "         MonteCarloArea\n"
    1909         "   ),\n"
    1910         "   ResultMC AS\n" // Group AnalysisMC by EnergySim Binning
    1911         "   (\n"
    1912         "      SELECT\n"
    1913         "         `.energysim`             AS `.energy`,\n"
    1914         "         ANY_VALUE(ThresholdW)    AS ThresholdW,\n"
    1915         "         ANY_VALUE(ThresholdW2)   AS ThresholdW2,\n"
    1916         "         ANY_VALUE(ThresholdN)    AS ThresholdN,\n"
    1917         "         ANY_VALUE(BiasSim)       AS Bias,\n"
    1918         "         ANY_VALUE(ResolutionSim) AS Resolution\n"
    1919         "      FROM\n"
    1920         "         AnalysisMC\n"
    1921         "      GROUP BY\n"
    1922         "         `.energy`\n"
    1923         "   )\n"
    1924         "   SELECT\n"
    1925         "      `.energy`, lo, hi,\n"
    1926         "      ThresholdW,\n"
    1927         "      SQRT(ThresholdW2) AS ErrThresholdW,\n"
    1928         "      ThresholdN,\n"
    1929         "      SQRT(ThresholdN)  AS ErrThresholdN,\n"         // Scale for Theta-Weights
    1930         "      ThresholdW        * 1000/(POW(10,hi)-POW(10,lo)) / Area / CountSim*ObsTime AS Flux,\n"
    1931         "      SQRT(ThresholdW2) * 1000/(POW(10,hi)-POW(10,lo)) / Area / CountSim*ObsTime AS ErrFlux,\n"
    1932         "      ThresholdN/CountN AS CutEfficiencyN,\n"
    1933         "      ThresholdW/CountW AS CutEfficiencyW,\n"
    1934         "      ThresholdN/CountN*(1/ThresholdN + 1/CountN) AS ErrCutEfficiencyN,\n"
    1935         "      ThresholdW/CountW*(ThresholdW2/POW(ThresholdW,2) + CountW2/POW(CountW,2)) AS ErrCutEfficiencyW,\n"
    1936         "      Bias,\n"
    1937         "      Resolution\n"
    1938         "   FROM\n"
    1939         "      ResultMC\n"
    1940         "   INNER JOIN\n"
    1941         "      Triggered USING (`.energy`)\n"
    1942         "   INNER JOIN\n"
    1943         "      BinningEnergySim ON `.energy`=bin\n"
    1944         "   CROSS JOIN\n"
    1945         "      ThetaSums, Area\n"
    1946         "   WHERE\n"
    1947         "      ThresholdW>0 AND ThresholdW2>0\n"
    1948         "   ORDER BY\n"
    1949         "      `.energy`\n"
    1950         ")";
    1951 
    1952     query9.parse();
    1953     for (auto it=env.cbegin(); it!=env.cend(); it++)
    1954         query9.template_defaults[it->first.c_str()] = it->second.c_str();
    1955 
    1956     //query9.template_defaults["area"] = area;
    1957     //query9.template_defaults["ontime"] = resX[0]["OnTime"].data();
    1958     //query9.template_defaults["count"]  = resX[0]["CountN"].data();
    1959 
    1960     if (print_queries)
    1961         PrintQuery(query9.str());
    1962 
    1963     qlog << query9 << ";\n" << endl;
    1964     if (connection.connected())
    1965     {
    1966         cout << query9.execute().info() << endl;
    1967         ShowWarnings(connection);
    1968         Dump(flog, connection, "Threshold");
    1969 
    1970         const auto sec9 = Time().UnixTime()-start9.UnixTime();
    1971         cout << "Execution time: " << sec9 << "s\n";
    1972 
    1973         // --------------------------------------------------------------------------
    1974 
    1975         const mysqlpp::StoreQueryResult res9 = connection.query("SELECT * FROM Threshold").store();
    1976 
    1977         sindent indentm(mlog, 4);
    1978 
    1979         mlog << '\n';
    1980         mlog << "// Energy Threshold\n";
    1981         mlog << "TGraphErrors g;\n";
    1982         mlog << "g.SetNameTitle(\"Threshold\", \"Energy Threshold\");\n";
    1983         for (auto ir=res9.cbegin(); ir!=res9.cend(); ir++)
    1984         {
    1985             // This is not efficient but easier to read and safer
    1986             const mysqlpp::Row &row = *ir;
    1987 
    1988             const double hi     = row["hi"];
    1989             const double lo     = row["lo"];
    1990             const double center = (hi+lo)/2;
    1991             const double width  = pow(10, hi)-pow(10, lo);
    1992             const double flux   = row["Flux"]   /width;
    1993             const double error  = row["ErrFlux"]/width;
    1994 
    1995             mlog << "g.SetPoint(g.GetN(), " << setw(8) << center << ", " << flux << ");\n";
    1996             mlog << "g.SetPointError(g.GetN()-1, 0, " << error << ");\n";
    1997         }
    1998         mlog << "g.GetXaxis()->SetTitle(\"lg(E/TeV)\");\n";
    1999         mlog << "g.GetYaxis()->SetTitle(\"dN/dE [cm^{-2} s^{-1} TeV^{-1}]\");\n";
    2000         mlog << indent(0) << "}" << endl;
    2001     }
    2002 
    2003     WriteHistogram(connection, {
    2004              .name     = "SimExcessW",
    2005              .title    = "Weighted Simulated Excess",
    2006              .dir      = "Esim/Simulation/Weighted",
    2007              .binningx = binning_esim,
    2008              .binningy = {},
    2009              .table    = "Threshold",
    2010              .x        = ".energy",
    2011              .y        = "",
    2012              .v        = "ThresholdW",
    2013              .err      = "ErrThresholdW",
    2014              .axisx    = "Energy lg(E_{sim}/GeV)",
    2015              .axisy    = "Weighted Counts",
    2016              .axisz    = "",
    2017              .stats    = true
    2018          });
    2019 
    2020     WriteHistogram(connection, {
    2021              .name     = "SimExcessN",
    2022              .title    = "Simulated Excess",
    2023              .dir      = "Esim/Simulation/Counts",
    2024              .binningx = binning_esim,
    2025              .binningy = {},
    2026              .table    = "Threshold",
    2027              .x        = ".energy",
    2028              .y        = "",
    2029              .v        = "ThresholdN",
    2030              .err      = "ErrThresholdN",
    2031              .axisx    = "Energy lg(E_{sim}/GeV)",
    2032              .axisy    = "Counts",
    2033              .axisz    = "",
    2034              .stats    = true
    2035          });
    2036 
    2037     WriteHistogram(connection, {
    2038              .name     = "SimSpectrumW",
    2039              .title    = "Weighted Simulated Excess Spectrum",
    2040              .dir      = "Esim/Simulation/Weighted",
    2041              .binningx = binning_esim,
    2042              .binningy = {},
    2043              .table    = "Threshold",
    2044              .x        = ".energy",
    2045              .y        = "",
    2046              .v        = "Flux",
    2047              .err      = "ErrFlux",
    2048              .axisx    = "Energy lg(E_{sim}/GeV)",
    2049              .axisy    = "dN/dE [cm^{-2} s^{-1} TeV^{-1}]",
    2050              .axisz    = "",
    2051              .stats    = true
    2052          });
    2053 
    2054     WriteHistogram(connection, {
    2055              .name     = "CutEfficiencyN",
    2056              .title    = "Cut Efficency (Unweighted)",
    2057              .dir      = "Esim/CutEfficiency",
    2058              .binningx = binning_esim,
    2059              .binningy = {},
    2060              .table    = "Threshold",
    2061              .x        = ".energy",
    2062              .y        = "",
    2063              .v        = "CutEfficiencyN",
    2064              .err      = "ErrCutEfficiencyN",
    2065              .axisx    = "Energy lg(E_{sim}/GeV)",
    2066              .axisy    = "Efficiency",
    2067              .axisz    = "",
    2068              .stats    = false
    2069          });
    2070 
    2071     WriteHistogram(connection, {
    2072              .name     = "CutEfficiencyW",
    2073              .title    = "Cut Efficency (Weighted)",
    2074              .dir      = "Esim/CutEfficiency",
    2075              .binningx = binning_esim,
    2076              .binningy = {},
    2077              .table    = "Threshold",
    2078              .x        = ".energy",
    2079              .y        = "",
    2080              .v        = "CutEfficiencyW",
    2081              .err      = "ErrCutEfficiencyW",
    2082              .axisx    = "Energy lg(E_{sim}/GeV)",
    2083              .axisy    = "Efficiency",
    2084              .axisz    = "",
    2085              .stats    = false
    2086          });
    2087 
    2088     WriteHistogram(connection, {
    2089              .name     = "Bias",
    2090              .title    = "Energy Bias",
    2091              .dir      = "Esim",
    2092              .binningx = binning_esim,
    2093              .binningy = {},
    2094              .table    = "Threshold",
    2095              .x        = ".energy",
    2096              .y        = "",
    2097              .v        = "Bias",
    2098              .err      = "Resolution",
    2099              .axisx    = "Energy lg(E_{sim}/GeV)",
    2100              .axisy    = "lg(E_{est}/E_{sim})",
    2101              .axisz    = "",
    2102              .stats    = false
    2103          });
    2104 
    2105     WriteHistogram(connection, {
    2106              .name     = "Resolution",
    2107              .title    = "Energy Resolution",
    2108              .dir      = "Esim",
    2109              .binningx = binning_esim,
    2110              .binningy = {},
    2111              .table    = "Threshold",
    2112              .x        = ".energy",
    2113              .y        = "",
    2114              .v        = "Resolution",
    2115              .err      = "",
    2116              .axisx    = "Energy lg(E_{sim}/GeV)",
    2117              .axisy    = "#sigma(lg(E_{est}/E_{sim}))",
    2118              .axisz    = "",
    2119              .stats    = false
    2120          });
    2121 
    2122 
    2123     // -------------------------------------------------------------------
    2124     // ---------------------------- Migration ----------------------------
    2125     // -------------------------------------------------------------------
    2126 
    2127     cout << separator("Migration") << '\n';
    2128 
    2129     Time start10;
    2130 
    2131     /* 02:obtain-migration-matrix.sql */
    2132 
    2133     // This query gets the analysis results versus Estimated Energy from the combined table
    2134     mysqlpp::Query query10(&connection);
    2135     query10 <<
    2136         "CREATE TEMPORARY TABLE Migration ENGINE=Memory AS\n"
    2137         "(\n"
    2138         "   SELECT\n"
    2139         "      `.energyest`,\n"
    2140         "      `.energysim`,\n"
    2141         "      BinningEnergySim.lo   AS EsimLo,\n"
    2142         "      BinningEnergySim.hi   AS EsimHi,\n"
    2143         "      BinningEnergyEst.lo   AS EestLo,\n"
    2144         "      BinningEnergyEst.hi   AS EestHi,\n"
    2145         "      ANY_VALUE(MigrationW) AS MigrationW,\n"
    2146         "      ANY_VALUE(MigrationN) AS MigrationN\n"
    2147         //     FIXME: Errors
    2148         "   FROM\n"
    2149         "      AnalysisMC\n"
    2150         "   INNER JOIN\n"
    2151         "      BinningEnergyEst ON `.energyest`=BinningEnergyEst.bin\n"
    2152         "   INNER JOIN\n"
    2153         "      BinningEnergySim ON `.energysim`=BinningEnergySim.bin\n"
    2154         "   GROUP BY\n"
    2155         "      `.energyest`, `.energysim`\n"
    2156         "   ORDER BY\n"
    2157         "      `.energyest`, `.energysim`\n"
    2158         ")";
    2159 
    2160     if (print_queries)
    2161         PrintQuery(query10.str());
    2162 
    2163     qlog << query10 << ";\n" << endl;
    2164     if (connection.connected())
    2165     {
    2166         cout << query10.execute().info() << endl;
    2167         ShowWarnings(connection);
    2168         Dump(flog, connection, "Migration");
    2169 
    2170         const auto sec10 = Time().UnixTime()-start10.UnixTime();
    2171         cout << "Execution time: " << sec10 << "s\n\n";
    2172     }
    2173 
    2174     WriteHistogram(connection, {
    2175              .name     = "MigrationN",
    2176              .title    = "Energy Migration",
    2177              .dir      = "",
    2178              .binningx = binning_esim,
    2179              .binningy = binning_eest,
    2180              .table    = "Migration",
    2181              .x        = ".energysim",
    2182              .y        = ".energyest",
    2183              .v        = "MigrationN",
    2184              .err      = "",
    2185              .axisx    = "Energy lg(E_{sim}/GeV)",
    2186              .axisy    = "Energy lg(E_{est}/GeV)",
    2187              .axisz    = "Counts",
    2188              .stats    = false
    2189          });
    2190 
    2191     WriteHistogram(connection, {
    2192              .name     = "MigrationW",
    2193              .title    = "Energy Migration",
    2194              .dir      = "",
    2195              .binningx = binning_esim,
    2196              .binningy = binning_eest,
    2197              .table    = "Migration",
    2198              .x        = ".energysim",
    2199              .y        = ".energyest",
    2200              .v        = "MigrationW",
    2201              .err      = "",
    2202              .axisx    = "Energy lg(E_{sim}/GeV)",
    2203              .axisy    = "Energy lg(E_{est}/GeV)",
    2204              .axisz    = "Weighted Counts",
    2205              .stats    = false
    2206          });
    2207 
    2208 
    2209 
    2210     // -------------------------------------------------------------------
    2211     // --------------------------- 11: Summary ---------------------------
     1993    }
     1994
     1995    // -------------------------------------------------------------------
     1996    // ----------------------------- Summary -----------------------------
     1997    // -------------------------------------------------------------------
    22121998
    22131999    cout << separator("Summary") << '\n';
Note: See TracChangeset for help on using the changeset viewer.