Changeset 19889


Ignore:
Timestamp:
12/12/19 16:25:49 (5 years ago)
Author:
tbretz
Message:
Added counting triggered events (Triggered) to display cut efficiencies, added total efficiencies, renamed ThetaHist to ThetaDistribution
File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/FACT++/src/spectrum.cc

    r19888 r19889  
    691691
    692692    // -------------------------------------------------------------------
    693     // --------------------------- 0: Binnings ---------------------------
     693    // ---------------------------- Binnings -----------------------------
     694    // -------------------------------------------------------------------
    694695
    695696    cout << separator("Binnings") << '\n';
    … …  
    704705
    705706    // -------------------------------------------------------------------
    706     // --------------------------- 1: DataFiles --------------------------
     707    // ---------------------------- DataFiles ----------------------------
     708    // -------------------------------------------------------------------
     709
    707710    cout << separator("DataFiles") << '\n';
    708711
    … …  
    754757
    755758    // -------------------------------------------------------------------
    756     // ------------------------ 2: ObservationTime -----------------------
     759    // ------------------------- ObservationTime -------------------------
     760    // -------------------------------------------------------------------
    757761
    758762    cout << separator("ObservationTime") << '\n';
    … …  
    867871
    868872    // -------------------------------------------------------------------
    869     // ------------------------ 3: MonteCarloFiles -----------------------
     873    // -------------------------- MonteCarloFiles ------------------------
     874    // -------------------------------------------------------------------
    870875
    871876    cout << separator("MonteCarloFiles") << '\n';
    … …  
    931936
    932937    // -------------------------------------------------------------------
    933     // ----------------------- 3b: Monte Carlo Area ----------------------
     938    // ------------------------- Monte Carlo Area ------------------------
     939    // -------------------------------------------------------------------
    934940
    935941    cout << separator("MonteCarloArea") << '\n';
    … …  
    982988
    983989    // -------------------------------------------------------------------
    984     // -------------------------- 4: ThetaHist ---------------------------
    985 
    986     cout << separator("ThetaHist") << '\n';
     990    // ------------------------ ThetaDistribution ------------------------
     991    // -------------------------------------------------------------------
     992
     993    cout << separator("ThetaDistribution") << '\n';
    987994
    988995    Time start4;
    … …  
    991998    mysqlpp::Query query4(&connection);
    992999    query4 <<
    993         "CREATE TEMPORARY TABLE ThetaHist\n"
     1000        "CREATE TEMPORARY TABLE ThetaDistribution\n"
    9941001        "(\n"
    9951002        "   `.theta`    SMALLINT UNSIGNED NOT NULL,\n"
    … …  
    10481055        cout << query4.execute().info() << endl;
    10491056        ShowWarnings(connection);
    1050         Dump(flog, connection, "ThetaHist");
     1057        Dump(flog, connection, "ThetaDistribution");
    10511058
    10521059        const auto sec4 = Time().UnixTime()-start4.UnixTime();
    … …  
    10561063        if (verbose>0)
    10571064        {
    1058             const mysqlpp::StoreQueryResult res4 = connection.query("SELECT * FROM ThetaHist").store();
     1065            const mysqlpp::StoreQueryResult res4 = connection.query("SELECT * FROM ThetaDistribution").store();
    10591066
    10601067            cout << "  Bin   MC Counts     OnTime ZdWeight\n";
    … …  
    10771084             .binningx = binning_theta,
    10781085             .binningy = {},
    1079              .table    = "ThetaHist",
     1086             .table    = "ThetaDistribution",
    10801087             .x        = ".theta",
    10811088             .y        = "",
    … …  
    10941101             .binningx = binning_theta,
    10951102             .binningy = {},
    1096              .table    = "ThetaHist",
     1103             .table    = "ThetaDistribution",
    10971104             .x        = ".theta",
    10981105             .y        = "",
    … …  
    11111118             .binningx = binning_theta,
    11121119             .binningy = {},
    1113              .table    = "ThetaHist",
     1120             .table    = "ThetaDistribution",
    11141121             .x        = ".theta",
    11151122             .y        = "",
    … …  
    11251132    // -------------------------------------------------------------------
    11261133    // ------------------------- 5: AnalysisData -------------------------
     1134    // -------------------------------------------------------------------
    11271135
    11281136    cout << separator("AnalysisData") << '\n';
    … …  
    12011209
    12021210    // -------------------------------------------------------------------
    1203     // --------------------------- 6: ResultMC ---------------------------
     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    // -------------------------------------------------------------------
    12041284
    12051285    cout << separator("ResultMC") << '\n';
    … …  
    12781358    }
    12791359
    1280 
    1281     // -------------------------------------------------------------------
    1282     // ----------------------- 7: SimulatedSpectrum ----------------------
     1360    // -------------------------------------------------------------------
     1361    // ------------------------- SimulatedSpectrum -----------------------
     1362    // -------------------------------------------------------------------
    12831363
    12841364    cout << separator("SimulatedSpectrum") << '\n';
    … …  
    13611441
    13621442    // -------------------------------------------------------------------
    1363     // --------------------------- 8: Spectrum ---------------------------
     1443    // ----------------------------- Spectrum ----------------------------
     1444    // -------------------------------------------------------------------
    13641445
    13651446    cout << separator("Spectrum") << '\n';
    … …  
    14091490        "         SUM(OnTime) AS ObsTime\n"
    14101491        "      FROM\n"
    1411         "         ThetaHist\n"
     1492        "         ThetaDistribution\n"
    14121493        "   ),\n"
    14131494        "   Area AS\n"
    … …  
    14621543        "      ( POW(ResultMC.ErrExcessW/ResultMC.ExcessW, 2) + POW(SQRT(SimulatedSpectrum.CountW2)/SimulatedSpectrum.CountW, 2) )\n"
    14631544        "         * POW(ResultMC.ExcessW/SimulatedSpectrum.CountW, 2) * CountSim/ObsTime AS ErrEfficiencyW,\n"
    1464         "      ( POW(ResultMC.ErrExcessN, 2) + POW(ResultMC.ExcessN, 2)/SimulatedSpectrum.CountN)/POW(SimulatedSpectrum.CountN, 2) AS ErrEfficiencyN\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"
    14651551        "   FROM\n"
    14661552        "      AnalysisData\n"
    … …  
    17401826
    17411827    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, {
    17421862             .name     = "Spectrum",
    17431863             .title    = "Differential Energy Spectrum",
    … …  
    17581878
    17591879    // -------------------------------------------------------------------
    1760     // -------------------------- 9: Threshold ---------------------------
     1880    // ---------------------------- Threshold ----------------------------
     1881    // -------------------------------------------------------------------
    17611882
    17621883    cout << separator("Threshold") << '\n';
    … …  
    17781899        "         SUM(OnTime) AS ObsTime\n"
    17791900        "      FROM\n"
    1780         "         ThetaHist\n"
     1901        "         ThetaDistribution\n"
    17811902        "   ),\n"
    17821903        "   Area AS\n"
    … …  
    18091930        "      ThresholdW        * 1000/(POW(10,hi)-POW(10,lo)) / Area / CountSim*ObsTime AS Flux,\n"
    18101931        "      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"
    18111936        "      Bias,\n"
    18121937        "      Resolution\n"
    18131938        "   FROM\n"
    18141939        "      ResultMC\n"
     1940        "   INNER JOIN\n"
     1941        "      Triggered USING (`.energy`)\n"
    18151942        "   INNER JOIN\n"
    18161943        "      BinningEnergySim ON `.energy`=bin\n"
    … …  
    19262053
    19272054    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, {
    19282089             .name     = "Bias",
    19292090             .title    = "Energy Bias",
    … …  
    19612122
    19622123    // -------------------------------------------------------------------
    1963     // -------------------------- 10: Migration --------------------------
     2124    // ---------------------------- Migration ----------------------------
     2125    // -------------------------------------------------------------------
    19642126
    19652127    cout << separator("Migration") << '\n';
Note: See TracChangeset for help on using the changeset viewer.