Index: trunk/FACT++/src/spectrum.cc =================================================================== --- trunk/FACT++/src/spectrum.cc (revision 19888) +++ trunk/FACT++/src/spectrum.cc (revision 19889) @@ -691,5 +691,6 @@ // ------------------------------------------------------------------- - // --------------------------- 0: Binnings --------------------------- + // ---------------------------- Binnings ----------------------------- + // ------------------------------------------------------------------- cout << separator("Binnings") << '\n'; @@ -704,5 +705,7 @@ // ------------------------------------------------------------------- - // --------------------------- 1: DataFiles -------------------------- + // ---------------------------- DataFiles ---------------------------- + // ------------------------------------------------------------------- + cout << separator("DataFiles") << '\n'; @@ -754,5 +757,6 @@ // ------------------------------------------------------------------- - // ------------------------ 2: ObservationTime ----------------------- + // ------------------------- ObservationTime ------------------------- + // ------------------------------------------------------------------- cout << separator("ObservationTime") << '\n'; @@ -867,5 +871,6 @@ // ------------------------------------------------------------------- - // ------------------------ 3: MonteCarloFiles ----------------------- + // -------------------------- MonteCarloFiles ------------------------ + // ------------------------------------------------------------------- cout << separator("MonteCarloFiles") << '\n'; @@ -931,5 +936,6 @@ // ------------------------------------------------------------------- - // ----------------------- 3b: Monte Carlo Area ---------------------- + // ------------------------- Monte Carlo Area ------------------------ + // ------------------------------------------------------------------- cout << separator("MonteCarloArea") << '\n'; @@ -982,7 +988,8 @@ // ------------------------------------------------------------------- - // -------------------------- 4: ThetaHist --------------------------- - - cout << separator("ThetaHist") << '\n'; + // ------------------------ ThetaDistribution ------------------------ + // ------------------------------------------------------------------- + + cout << separator("ThetaDistribution") << '\n'; Time start4; @@ -991,5 +998,5 @@ mysqlpp::Query query4(&connection); query4 << - "CREATE TEMPORARY TABLE ThetaHist\n" + "CREATE TEMPORARY TABLE ThetaDistribution\n" "(\n" " `.theta` SMALLINT UNSIGNED NOT NULL,\n" @@ -1048,5 +1055,5 @@ cout << query4.execute().info() << endl; ShowWarnings(connection); - Dump(flog, connection, "ThetaHist"); + Dump(flog, connection, "ThetaDistribution"); const auto sec4 = Time().UnixTime()-start4.UnixTime(); @@ -1056,5 +1063,5 @@ if (verbose>0) { - const mysqlpp::StoreQueryResult res4 = connection.query("SELECT * FROM ThetaHist").store(); + const mysqlpp::StoreQueryResult res4 = connection.query("SELECT * FROM ThetaDistribution").store(); cout << " Bin MC Counts OnTime ZdWeight\n"; @@ -1077,5 +1084,5 @@ .binningx = binning_theta, .binningy = {}, - .table = "ThetaHist", + .table = "ThetaDistribution", .x = ".theta", .y = "", @@ -1094,5 +1101,5 @@ .binningx = binning_theta, .binningy = {}, - .table = "ThetaHist", + .table = "ThetaDistribution", .x = ".theta", .y = "", @@ -1111,5 +1118,5 @@ .binningx = binning_theta, .binningy = {}, - .table = "ThetaHist", + .table = "ThetaDistribution", .x = ".theta", .y = "", @@ -1125,4 +1132,5 @@ // ------------------------------------------------------------------- // ------------------------- 5: AnalysisData ------------------------- + // ------------------------------------------------------------------- cout << separator("AnalysisData") << '\n'; @@ -1201,5 +1209,77 @@ // ------------------------------------------------------------------- - // --------------------------- 6: ResultMC --------------------------- + // ----------------------------- Triggered --------------------------- + // ------------------------------------------------------------------- + + cout << separator("Triggered") << '\n'; + + Time start6a; + + /* 02:analyze-simulation.sql */ + + // This table combines the analysis results vs. Binning in Estimated Energy and Simulated Energy + mysqlpp::Query query6a(&connection); + query6a << + "CREATE TEMPORARY TABLE Triggered\n" + "(\n" + " `.energy` SMALLINT UNSIGNED NOT NULL,\n" + " CountN DOUBLE NOT NULL,\n" + " CountW DOUBLE NOT NULL,\n" + " CountW2 DOUBLE NOT NULL,\n" + " INDEX (`.energy`) USING HASH\n" + ") ENGINE=Memory\n" + "AS\n" + "(\n" + " WITH TriggeredEvents AS\n" + " (\n" + " SELECT\n" + " INTERVAL(Zd, %100:theta) AS `.theta`,\n" + " INTERVAL(log10(Energy), %101:energy) AS `.energy`,\n" + " (%102:spectrum)/pow(Energy, SpectralIndex) AS SpectralWeight\n" + " FROM\n" + " MonteCarloFiles\n" + " LEFT JOIN\n" + " factmc.RunInfoMC USING (FileId)\n" + " LEFT JOIN\n" + " factmc.EventsMC USING (FileId)\n" + " )\n" + " SELECT\n" + " `.energy`,\n" + " COUNT(*) AS CountN,\n" + " SUM(ZdWeight*SpectralWeight) AS CountW,\n" + " SUM(POW(ErrZdWeight*SpectralWeight,2)) AS `CountW2`\n" + " FROM\n" + " TriggeredEvents\n" + " INNER JOIN\n" + " ThetaDistribution USING(`.theta`)\n" + " GROUP BY\n" + " `.energy`\n" + ")"; + + query6a.parse(); + for (auto it=env.cbegin(); it!=env.cend(); it++) + query6a.template_defaults[it->first.c_str()] = it->second.c_str(); + + query6a.template_defaults["energy"] = strth.c_str(); + query6a.template_defaults["theta"] = strzd.c_str(); + query6a.template_defaults["spectrum"] = spectrum.c_str(); + + if (print_queries) + PrintQuery(query6a.str()); + + qlog << query6a << ";\n" << endl; + if (connection.connected()) + { + cout << query6a.execute().info() << endl; + ShowWarnings(connection); + Dump(flog, connection, "Triggered"); + + const auto sec6a = Time().UnixTime()-start6a.UnixTime(); + cout << "Execution time: " << sec6a << "s\n\n"; + } + + // ------------------------------------------------------------------- + // ----------------------------- ResultMC ---------------------------- + // ------------------------------------------------------------------- cout << separator("ResultMC") << '\n'; @@ -1278,7 +1358,7 @@ } - - // ------------------------------------------------------------------- - // ----------------------- 7: SimulatedSpectrum ---------------------- + // ------------------------------------------------------------------- + // ------------------------- SimulatedSpectrum ----------------------- + // ------------------------------------------------------------------- cout << separator("SimulatedSpectrum") << '\n'; @@ -1361,5 +1441,6 @@ // ------------------------------------------------------------------- - // --------------------------- 8: Spectrum --------------------------- + // ----------------------------- Spectrum ---------------------------- + // ------------------------------------------------------------------- cout << separator("Spectrum") << '\n'; @@ -1409,5 +1490,5 @@ " SUM(OnTime) AS ObsTime\n" " FROM\n" - " ThetaHist\n" + " ThetaDistribution\n" " ),\n" " Area AS\n" @@ -1462,5 +1543,10 @@ " ( POW(ResultMC.ErrExcessW/ResultMC.ExcessW, 2) + POW(SQRT(SimulatedSpectrum.CountW2)/SimulatedSpectrum.CountW, 2) )\n" " * POW(ResultMC.ExcessW/SimulatedSpectrum.CountW, 2) * CountSim/ObsTime AS ErrEfficiencyW,\n" - " ( POW(ResultMC.ErrExcessN, 2) + POW(ResultMC.ExcessN, 2)/SimulatedSpectrum.CountN)/POW(SimulatedSpectrum.CountN, 2) AS ErrEfficiencyN\n" + " ( POW(ResultMC.ErrExcessN, 2) + POW(ResultMC.ExcessN, 2)/SimulatedSpectrum.CountN)/POW(SimulatedSpectrum.CountN, 2) AS ErrEfficiencyN,\n" + " Area/10000*ResultMC.ExcessW/SimulatedSpectrum.CountW * CountSim/ObsTime AS EffAreaW,\n" + " Area/10000*ResultMC.ExcessN/SimulatedSpectrum.CountN AS EffAreaN,\n" + " Area/10000*( POW(ResultMC.ErrExcessW/ResultMC.ExcessW, 2) + POW(SQRT(SimulatedSpectrum.CountW2)/SimulatedSpectrum.CountW, 2) )\n" + " * POW(ResultMC.ExcessW/SimulatedSpectrum.CountW, 2) * CountSim/ObsTime AS ErrEffAreaW,\n" + " Area/10000*( POW(ResultMC.ErrExcessN, 2) + POW(ResultMC.ExcessN, 2)/SimulatedSpectrum.CountN)/POW(SimulatedSpectrum.CountN, 2) AS ErrEffAreaN\n" " FROM\n" " AnalysisData\n" @@ -1740,4 +1826,38 @@ WriteHistogram(connection, { + .name = "EffAreaN", + .title = "Effective Area (Counts)", + .dir = "Eest/EffArea", + .binningx = binning_eest, + .binningy = {}, + .table = "Spectrum", + .x = ".energy", + .y = "", + .v = "EffAreaN", + .err = "ErrEffAreaN", + .axisx = "Energy lg(E_{sim}/GeV)", + .axisy = "Effective Area A [m^{2}]", + .axisz = "", + .stats = false + }); + + WriteHistogram(connection, { + .name = "EffAreaW", + .title = "Effective Area (Weighted)", + .dir = "Eest/EffArea", + .binningx = binning_eest, + .binningy = {}, + .table = "Spectrum", + .x = ".energy", + .y = "", + .v = "EffAreaW", + .err = "ErrEffAreaW", + .axisx = "Energy lg(E_{sim}/GeV)", + .axisy = "Effective Area A [m^{2}]", + .axisz = "", + .stats = false + }); + + WriteHistogram(connection, { .name = "Spectrum", .title = "Differential Energy Spectrum", @@ -1758,5 +1878,6 @@ // ------------------------------------------------------------------- - // -------------------------- 9: Threshold --------------------------- + // ---------------------------- Threshold ---------------------------- + // ------------------------------------------------------------------- cout << separator("Threshold") << '\n'; @@ -1778,5 +1899,5 @@ " SUM(OnTime) AS ObsTime\n" " FROM\n" - " ThetaHist\n" + " ThetaDistribution\n" " ),\n" " Area AS\n" @@ -1809,8 +1930,14 @@ " ThresholdW * 1000/(POW(10,hi)-POW(10,lo)) / Area / CountSim*ObsTime AS Flux,\n" " SQRT(ThresholdW2) * 1000/(POW(10,hi)-POW(10,lo)) / Area / CountSim*ObsTime AS ErrFlux,\n" + " ThresholdN/CountN AS CutEfficiencyN,\n" + " ThresholdW/CountW AS CutEfficiencyW,\n" + " ThresholdN/CountN*(1/ThresholdN + 1/CountN) AS ErrCutEfficiencyN,\n" + " ThresholdW/CountW*(ThresholdW2/POW(ThresholdW,2) + CountW2/POW(CountW,2)) AS ErrCutEfficiencyW,\n" " Bias,\n" " Resolution\n" " FROM\n" " ResultMC\n" + " INNER JOIN\n" + " Triggered USING (`.energy`)\n" " INNER JOIN\n" " BinningEnergySim ON `.energy`=bin\n" @@ -1926,4 +2053,38 @@ WriteHistogram(connection, { + .name = "CutEfficiencyN", + .title = "Cut Efficency (Unweighted)", + .dir = "Esim/CutEfficiency", + .binningx = binning_esim, + .binningy = {}, + .table = "Threshold", + .x = ".energy", + .y = "", + .v = "CutEfficiencyN", + .err = "ErrCutEfficiencyN", + .axisx = "Energy lg(E_{sim}/GeV)", + .axisy = "Efficiency", + .axisz = "", + .stats = false + }); + + WriteHistogram(connection, { + .name = "CutEfficiencyW", + .title = "Cut Efficency (Weighted)", + .dir = "Esim/CutEfficiency", + .binningx = binning_esim, + .binningy = {}, + .table = "Threshold", + .x = ".energy", + .y = "", + .v = "CutEfficiencyW", + .err = "ErrCutEfficiencyW", + .axisx = "Energy lg(E_{sim}/GeV)", + .axisy = "Efficiency", + .axisz = "", + .stats = false + }); + + WriteHistogram(connection, { .name = "Bias", .title = "Energy Bias", @@ -1961,5 +2122,6 @@ // ------------------------------------------------------------------- - // -------------------------- 10: Migration -------------------------- + // ---------------------------- Migration ---------------------------- + // ------------------------------------------------------------------- cout << separator("Migration") << '\n';