Index: trunk/FACT++/src/spectrum.cc =================================================================== --- trunk/FACT++/src/spectrum.cc (revision 19889) +++ trunk/FACT++/src/spectrum.cc (revision 19895) @@ -17,4 +17,6 @@ #include "TH1.h" #include "TH2.h" +#include "TRolke.h" +#include "TFeldmanCousins.h" #endif @@ -130,16 +132,17 @@ #endif , "Database link as in\n\tuser:password@server[:port]/database[?compress=0|1].") - ("out,o", var(conf.GetName()), "Defines the prefix (with path) of the output files.") - ("dry-run", po_switch(), "Only write the queries to the query.log file. Internally overwrites uri with an empty string.") + ("out,o", var(conf.GetName()), "Defines the prefix (with path) of the output files.") + ("confidence-level,c", var(0.99), "Confidence level for the calculation of the upper limits.") + ("feldman-cousins", po_bool(), "Calculate Feldman-Cousins ULs (slow and only minor difference to Rolke).") ; po::options_description binnings("Binnings"); binnings.add_options() - ("theta", var(Binning(90, 0, 90)), "Add equidistant bins in theta (degrees). Syntax: N,lo,hi (Number N of equidistant bins between lo and hi)") - ("theta-bin", vars(), "Add a bin-edge to the theta binning (degree)") - ("esim", var(Binning(15, 2, 5)), "Add equidistant bins in log10 simulated energy. Syntax: N,lo,hi (Number N of equidistant bins between lo and hi)") - ("esim-bin", vars(), "Add a bin-edge to the binnnig in log10 simulated enegry") - ("eest", var(Binning(15, 2, 5)), "Add equidistant bins in log10 estimated energy. Syntax: N,lo,hi (Number N of equidistant bins between lo and hi)") - ("eest-bin", vars(), "Add a bin-edge to the binning in log10 estimated enegry") + ("theta", var(Binning(90, 0, 90)), "Add equidistant bins in theta (degrees). Syntax: N,lo,hi (Number N of equidistant bins between lo and hi)") + ("theta-bin", vars(), "Add a bin-edge to the theta binning (degree)") + ("energy-dense", var(Binning(30, 2, 5)), "Add equidistant bins in log10 simulated energy. Syntax: N,lo,hi (Number N of equidistant bins between lo and hi)") + ("energy-dense-bin", vars(), "Add a bin-edge to the binnnig in log10 simulated enegry") + ("energy-sparse", var(Binning(15, 2, 5)), "Add equidistant bins in log10 estimated energy. Syntax: N,lo,hi (Number N of equidistant bins between lo and hi)") + ("energy-sparse-bin", vars(), "Add a bin-edge to the binning in log10 estimated enegry") ; @@ -147,6 +150,6 @@ queries.add_options() ("analysis", var("analysis.sql"), "File with the analysis query. A default file is created automatically in the directory it does not exist.") - ("data", var("data.sql"), "File with the query which creates the data summary. A default file is created automatically in the directory it does not exist.") - ("simulation", var("simulation.sql"), "File with the query which creates the simulation summary. A default file is created automatically in the directory it does not exist.") + //("data", var("data.sql"), "File with the query which creates the data summary. A default file is created automatically in the directory it does not exist.") + //("simulation", var("simulation.sql"), "File with the query which creates the simulation summary. A default file is created automatically in the directory it does not exist.") ; @@ -162,10 +165,12 @@ po::options_description debug("Debug options"); debug.add_options() - ("print-connection", po_switch(), "Print database connection information") + ("dry-run", po_bool(), "Only write the queries to the query.log file. Internally overwrites uri with an empty string.") + ("print-connection", po_bool(), "Print database connection information") #ifdef HAVE_HIGHLIGHT - ("print-queries", po_switch(), "Print all queries to the console. They are automatically highlighted.") + ("print-queries", po_bool(), "Print all queries to the console. They are automatically highlighted.") #else - ("print-queries", po_switch(), "Print all queries to the console. (For highlighting recompile with 'libsource-highlight-dev' installed)") + ("print-queries", po_bool(), "Print all queries to the console. (For highlighting recompile with 'libsource-highlight-dev' installed)") #endif + ("mc-only", po_bool(), "Do not run an data realated queries (except observation times)") ("verbose,v", var(1), "Verbosity (0: quiet, 1: default, 2: more, 3, ...)") ; @@ -327,4 +332,7 @@ srchilite::SourceHighlight sourceHighlight("esc256.outlang"); sourceHighlight.setStyleFile("esc256.style"); + sourceHighlight.setGenerateLineNumbers(); + sourceHighlight.setLineNumberDigits(3); + //sourceHighlight.setLineNumberPad(' ') sourceHighlight.highlight(qstr, cout, "sql.lang"); cout << endl; @@ -341,6 +349,6 @@ "(\n" " bin SMALLINT UNSIGNED NOT NULL,\n" - " lo FLOAT NOT NULL,\n" - " hi FLOAT NOT NULL,\n" + " lo DOUBLE NOT NULL,\n" + " hi DOUBLE NOT NULL,\n" " PRIMARY KEY (bin) USING HASH\n" ")"; @@ -349,6 +357,4 @@ if (connection.connected()) query0.execute(); - - //connection.query("ALTER TABLE Binning"+name+" AUTO_INCREMENT=0").execute(); mysqlpp::Query query1(&connection); @@ -426,4 +432,65 @@ }; +#ifdef HAVE_ROOT + +TH1 *CreateHistogram(const Histogram &hist) +{ + const char *name = hist.name.empty() ? hist.v.c_str() : hist.name.c_str(); + + cout << "Creating Histogram '" << hist.dir << "/" << name << "'" << endl; + + const auto vecx = hist.binningx.vec(); + const auto vecy = hist.binningy.vec(); + + TH1 *h = 0; + + if (hist.y.empty()) + { + h = hist.binningx.equidist ? + new TH1D(name, hist.title.c_str(), vecx.size()-1, vecx.front(), vecx.back()) : + new TH1D(name, hist.title.c_str(), vecx.size()-1, vecx.data()); + } + else + { + if (hist.binningy.equidist) + { + h = hist.binningx.equidist ? + new TH2D(name, hist.title.c_str(), vecx.size()-1, vecx.front(), vecx.back(), vecy.size()-1, vecy.front(), vecy.back()) : + new TH2D(name, hist.title.c_str(), vecx.size()-1, vecx.data(), vecy.size()-1, vecy.front(), vecy.back()); + } + else + { + h = hist.binningx.equidist ? + new TH2D(name, hist.title.c_str(), vecx.size()-1, vecx.front(), vecx.back(), vecy.size()-1, vecy.front(), vecy.back()) : + new TH2D(name, hist.title.c_str(), vecx.size()-1, vecx.data(), vecy.size()-1, vecy.front(), vecy.back()); + } + } + + h->SetXTitle(hist.axisx.c_str()); + h->SetYTitle(hist.axisy.c_str()); + h->SetZTitle(hist.axisz.c_str()); + + h->SetMarkerStyle(kFullDotMedium); + + h->SetStats(hist.stats); + + return h; +} + +void WriteHistogram(TH1 *h, const string &directory) +{ + gFile->cd(); + TDirectory *dir = gFile->GetDirectory(directory.c_str()); + if (dir) + dir->cd(); + else + { + gFile->mkdir(directory.c_str()); + gFile->cd(directory.c_str()); + } + h->Write(); +} +#endif + void WriteHistogram(Database &connection, const Histogram &hist) { @@ -431,4 +498,6 @@ if (!connection.connected()) return; + + TH1 *h = CreateHistogram(hist); mysqlpp::Query query(&connection); @@ -449,80 +518,56 @@ const mysqlpp::StoreQueryResult res = query.store(); - const auto vecx = hist.binningx.vec(); - const auto vecy = hist.binningy.vec(); - - TH1 *h = 0; - - if (hist.y.empty()) - { - h = hist.binningx.equidist ? - new TH1D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.front(), vecx.back()) : - new TH1D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.data()); - } - else - { - if (hist.binningy.equidist) - { - h = hist.binningx.equidist ? - new TH2D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.front(), vecx.back(), vecy.size()-1, vecy.front(), vecy.back()) : - new TH2D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.data(), vecy.size()-1, vecy.front(), vecy.back()); - } - else - { - h = hist.binningx.equidist ? - new TH2D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.front(), vecx.back(), vecy.size()-1, vecy.front(), vecy.back()) : - new TH2D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.data(), vecy.size()-1, vecy.front(), vecy.back()); - } - } - - h->SetXTitle(hist.axisx.c_str()); - h->SetYTitle(hist.axisy.c_str()); - h->SetZTitle(hist.axisz.c_str()); - - h->SetMarkerStyle(kFullDotMedium); - - h->SetStats(hist.stats); - for (auto ir=res.cbegin(); ir!=res.cend(); ir++) { const auto &row = *ir; - - const uint32_t x = row["X"]; - const double v = row["V"]; - if (x>uint32_t(h->GetNbinsX()+1)) - continue; try { - const uint32_t y = row["Y"]; - if (y>uint32_t(h->GetNbinsY()+1)) + const uint32_t x = row["X"]; + const double v = row["V"]; + if (x>uint32_t(h->GetNbinsX()+1)) continue; - h->SetBinContent(x, y, v); - - } - catch (const mysqlpp::BadFieldName &) - { - h->SetBinContent(x, v); try { - h->SetBinError(x, double(row["E"])); + const uint32_t y = row["Y"]; + if (y>uint32_t(h->GetNbinsY()+1)) + continue; + + h->SetBinContent(x, y, v); + } catch (const mysqlpp::BadFieldName &) { + h->SetBinContent(x, v); + try + { + h->SetBinError(x, double(row["E"])); + } + catch (const mysqlpp::BadFieldName &) + { + } } } - } - - gFile->cd(); - TDirectory *dir = gFile->GetDirectory(hist.dir.c_str()); - if (dir) - dir->cd(); - else - { - gFile->mkdir(hist.dir.c_str()); - gFile->cd(hist.dir.c_str()); - } - h->Write(); + catch (const mysqlpp::BadConversion &b) + { + cerr << b.what() << endl; + } + } + + WriteHistogram(h, hist.dir); + delete h; +#endif +} + +void WriteHistogram(Database &connection, const Histogram &hist, const map &data) +{ +#ifdef HAVE_ROOT + TH1 *h = CreateHistogram(hist); + + for (auto ir=data.cbegin(); ir!=data.cend(); ir++) + h->SetBinContent(ir->first, ir->second); + + WriteHistogram(h, hist.dir); delete h; #endif @@ -581,19 +626,22 @@ // ----------------------------- Evaluate options -------------------------- - const string uri = conf.Get("dry-run") ? "" : conf.Get("uri"); - const string out = conf.Get("out"); - const uint16_t verbose = conf.Get("verbose"); + const string uri = conf.Get("dry-run") ? "" : conf.Get("uri"); + const string out = conf.Get("out"); + const uint16_t verbose = conf.Get("verbose"); + const double confidence = conf.Get("confidence-level"); + const bool feldman = conf.Get("feldman-cousins"); const bool print_connection = conf.Get("print-connection"); const bool print_queries = conf.Get("print-queries"); - - Binning binning_theta = conf.Get("theta")+conf.Vec("theta-bin"); - Binning binning_esim = conf.Get("esim"); - Binning binning_eest = conf.Get("eest"); + const bool mc_only = conf.Get("mc-only"); + + Binning binning_theta = conf.Get("theta")+conf.Vec("theta-bin"); + Binning binning_dense = conf.Get("energy-dense"); + Binning binning_sparse = conf.Get("energy-sparse"); cout << '\n'; - cout << "Binning 'theta': " << binning_theta.str() << endl; - cout << "Binning 'Esim': " << binning_esim.str() << endl; - cout << "Binning 'Eest': " << binning_eest.str() << endl; + cout << "Binning 'theta': " << binning_theta.str() << endl; + cout << "Binning 'dense': " << binning_dense.str() << endl; + cout << "Binning 'sparse': " << binning_sparse.str() << endl; const uint16_t source_key = conf.Get("source-key"); @@ -604,12 +652,15 @@ cout << "\n"; - const string analysis_sql = CreateResource(conf, "analysis", "analysis.sql", RESOURCE(std::string, spectrum_analysis_sql)); - const string data_sql = CreateResource(conf, "data", "data.sql", RESOURCE(std::string, spectrum_data_sql)); - const string simulation_sql = CreateResource(conf, "simulation", "simulation.sql", RESOURCE(std::string, spectrum_simulation_sql)); + const string analysis_sql = CreateResource(conf, "analysis", "analysis.sql", RESOURCE(std::string, spectrum_analysis_sql)); + const string data_sql = RESOURCE(std::string, spectrum_data_sql); + const string simulation_sql = RESOURCE(std::string, spectrum_simulation_sql); + const string spectrum_sql = RESOURCE(std::string, spectrum_spectrum_sql); + const string summary_sim_sql = RESOURCE(std::string, spectrum_summary_sim_sql); + const string summary_est_sql = RESOURCE(std::string, spectrum_summary_est_sql); cout << endl; - const string strzd = binning_theta.list(); - const string stre = binning_eest.list(); - const string strth = binning_esim.list(); + const string str_theta = binning_theta.list(); + const string str_dense = binning_dense.list(); + const string str_sparse = binning_sparse.list(); // ------------------------------------------------------------------------- @@ -696,11 +747,11 @@ cout << separator("Binnings") << '\n'; - CreateBinning(connection, qlog, binning_theta, "Theta"); - CreateBinning(connection, qlog, binning_eest, "EnergyEst"); - CreateBinning(connection, qlog, binning_esim, "EnergySim"); + CreateBinning(connection, qlog, binning_theta, "Theta"); + CreateBinning(connection, qlog, binning_dense, "Energy_dense"); + CreateBinning(connection, qlog, binning_sparse, "Energy_sparse"); Dump(flog, connection, "BinningTheta"); - Dump(flog, connection, "BinningEnergyEst"); - Dump(flog, connection, "BinningEnergySim"); + Dump(flog, connection, "BinningEnergy_dense"); + Dump(flog, connection, "BinningEnergy_sparse"); // ------------------------------------------------------------------- @@ -719,5 +770,5 @@ " FileId INT UNSIGNED NOT NULL,\n" " PRIMARY KEY (FileId) USING HASH\n" - ") ENGINE=Memory\n" + ") ENGINE=MEMORY\n" "AS\n" "(\n" @@ -770,7 +821,7 @@ "(\n" " `.theta` SMALLINT UNSIGNED NOT NULL,\n" - " OnTime FLOAT NOT NULL,\n" + " OnTime DOUBLE NOT NULL,\n" " PRIMARY KEY (`.theta`) USING HASH\n" - ") ENGINE=Memory\n" + ") ENGINE=MEMORY\n" "AS\n" "(\n" @@ -792,5 +843,5 @@ query2.template_defaults[it->first.c_str()] = it->second.c_str(); - query2.template_defaults["bins"] = strzd.c_str(); + query2.template_defaults["bins"] = str_theta.c_str(); if (print_queries) @@ -808,66 +859,4 @@ } - /* - - SELECT - min(first_number) as first_number, - max(last_number) as last_number - from - ( - select - first_number, - last_number - sum(grp) over(order by first_number) as grp - from - ( - select - bin AS first_number, - LEAD(bin) AS last_number - IF(first_number <> lag(bin, 2) over (order by bin) + 1, 1, 0) as grp - from - t1 - ) - ) - group by grp - order by 1 - */ - - /* -WITH Hallo AS (WITH MyTest AS (WITH MyTable AS (SELECT bin FROM Test -WHERE cnt>0 GROUP BY bin ORDER BY bin) SELECT bin AS first_number, -LAG(bin) OVER(ORDER BY bin) AS prev, LEAD(bin) OVER(ORDER BY bin) AS -last_number, IF(bin <> (LAG(bin) OVER(ORDER BY bin)+1), 1, 0) AS -grp_prev, IF(bin <> (LEAD(bin) OVER(ORDER BY bin)-1), 1, 0) as grp_nxt -FROM MyTable ) SELECT first_number, last_number, sum(grp_nxt) -over(order by first_number) as grp FROM MyTest) SELECT * FROM Hallo - - - WITH MyTable AS -(SELECT bin FROM Test WHERE cnt>0 GROUP BY bin ORDER BY bin) -SELECT bin, LAG(bin) OVER(ORDER BY bin) AS prev, LEAD(bin) OVER(ORDER BY bin) AS nxt, -IF(bin = (LAG(bin) OVER(ORDER BY bin)+1), 1, 0) AS grp_prev, IF(bin = (LEAD(bin) OVER(ORDER BY bin)-1), 1, 0) as grp_nxt FROM MyTable - -"-1","0" -"0","0" -"2","3" 2-4 -"3","5" -"4","6" -"5","0" -"6","0" -"7","7" 7-7 -"8","0" -"9","9" 9-12 -"10","10" -"11","11" -"12","12" -"13","0" -"14","0" -"15","0" -"16","16" 16-19 -"17","17" -"18","18" -"19","19" - */ - // ------------------------------------------------------------------- // -------------------------- MonteCarloFiles ------------------------ @@ -878,5 +867,4 @@ Time start3; - /* 02:get-monte-carlo-files.sql */ mysqlpp::Query query3(&connection); query3 << @@ -885,5 +873,5 @@ " FileId INT UNSIGNED NOT NULL,\n" " PRIMARY KEY (FileId) USING HASH\n" - ") ENGINE=Memory\n" + ") ENGINE=MEMORY\n" "AS\n" "(\n" @@ -897,5 +885,4 @@ " factmc.RunInfoMC\n" " ON\n" - //" ThetaMin BETWEEN lo AND hi OR ThetaMax BETWEEN lo AND hi\n" // Includes BOTH edges " (ThetaMin>=lo AND ThetaMinlo AND ThetaMax<=hi)\n" " WHERE\n" @@ -905,15 +892,4 @@ " FileId\n" // In order: faster ")"; -/* - " SELECT\n" - " FileId\n" - " FROM\n" - " factmc.RunInfoMC\n" - " WHERE\n" - " PartId=1 AND\n" - " FileId%%2=0 AND\n" - " (%0:range)\n" - ")"; -*/ query3.parse(); @@ -941,10 +917,9 @@ cout << separator("MonteCarloArea") << '\n'; - Time start3b; - - /* 02:get-monte-carlo-files.sql */ - mysqlpp::Query query3b(&connection); - query3b << - "CREATE TEMPORARY TABLE MonteCarloArea ENGINE=Memory\n" + Time start4; + + mysqlpp::Query query4(&connection); + query4 << + "CREATE TEMPORARY TABLE MonteCarloArea ENGINE=MEMORY\n" "AS\n" "(\n" @@ -964,15 +939,15 @@ ")"; - query3b.parse(); + query4.parse(); for (auto it=env.cbegin(); it!=env.cend(); it++) - query3b.template_defaults[it->first.c_str()] = it->second.c_str(); + query4.template_defaults[it->first.c_str()] = it->second.c_str(); if (print_queries) - PrintQuery(query3b.str()); - - qlog << query3b << ";\n" << endl; + PrintQuery(query4.str()); + + qlog << query4 << ";\n" << endl; if (connection.connected()) { - cout << query3b.execute().info() << endl; + cout << query4.execute().info() << endl; ShowWarnings(connection); if (Dump(flog, connection, "MonteCarloArea")!=1) @@ -982,48 +957,134 @@ } - const auto sec3b = Time().UnixTime()-start3b.UnixTime(); - cout << "Execution time: " << sec3b << "s\n\n"; - } - - - // ------------------------------------------------------------------- - // ------------------------ ThetaDistribution ------------------------ - // ------------------------------------------------------------------- - - cout << separator("ThetaDistribution") << '\n'; - - Time start4; - - /* 02:get-theta-distribution.sql */ - mysqlpp::Query query4(&connection); - query4 << - "CREATE TEMPORARY TABLE ThetaDistribution\n" + const auto sec4 = Time().UnixTime()-start4.UnixTime(); + cout << "Execution time: " << sec4 << "s\n\n"; + } + + + // ------------------------------------------------------------------- + // ----------------------------- SummaryMC --------------------------- + // ------------------------------------------------------------------- + + cout << separator("SummaryOriginalMC") << '\n'; + + Time start5; + + // This table combines the analysis results vs. Binning in Estimated Energy and Simulated Energy + mysqlpp::Query query5(&connection); + query5 << + "CREATE TEMPORARY TABLE Summary%100:table\n" "(\n" - " `.theta` SMALLINT UNSIGNED NOT NULL,\n" - " lo DOUBLE NOT NULL COMMENT 'Lower edge of zenith distance bin in degree',\n" - " hi DOUBLE NOT NULL COMMENT 'Upper edge of zenith distance bin in degree',\n" - " CountN INT UNSIGNED NOT NULL,\n" - " ErrCountN DOUBLE NOT NULL,\n" - " OnTime FLOAT NOT NULL,\n" - " ZdWeight DOUBLE NOT NULL COMMENT 'tau(delta theta)',\n" - " ErrZdWeight DOUBLE NOT NULL COMMENT 'sigma(tau)',\n" - " PRIMARY KEY (`.theta`) USING HASH\n" - ") ENGINE=Memory\n" + " `.theta` SMALLINT UNSIGNED NOT NULL,\n" + " `.sparse_sim` SMALLINT UNSIGNED NOT NULL,\n" + " `.dense_sim` SMALLINT UNSIGNED NOT NULL,\n" + " CountN DOUBLE NOT NULL,\n" + " SumW DOUBLE NOT NULL,\n" + " SumW2 DOUBLE NOT NULL,\n" + " INDEX (`.theta`) USING HASH,\n" + " INDEX (`.sparse_sim`) USING HASH,\n" + " INDEX (`.dense_sim`) USING HASH\n" + ") ENGINE=MEMORY\n" "AS\n" "(\n" - " WITH EventCount AS\n" + " WITH Table0 AS\n" " (\n" " SELECT\n" - " INTERVAL(DEGREES(Theta), %100:bins) AS `.theta`,\n" - " COUNT(*) AS CountN\n" + " INTERVAL(%101:column, %102:theta) AS `.theta`,\n" + " INTERVAL(LOG10(Energy), %103:sparse) AS `.sparse_sim`,\n" + " INTERVAL(LOG10(Energy), %104:dense) AS `.dense_sim`,\n" + " (%105:spectrum)/pow(Energy, SpectralIndex) AS SpectralWeight\n" " FROM\n" " MonteCarloFiles\n" " LEFT JOIN\n" - " factmc.OriginalMC USING(FileId)\n" + " factmc.RunInfoMC USING (FileId)\n" + " LEFT JOIN\n" + " factmc.%100:table USING (FileId)\n" + " )\n" + " SELECT\n" + " `.theta`,\n" + " `.sparse_sim`,\n" + " `.dense_sim`,\n" + " COUNT(*) AS CountN,\n" + " SUM( SpectralWeight ) AS SumW,\n" + " SUM(POW(SpectralWeight,2)) AS SumW2\n" + " FROM\n" + " Table0\n" + " GROUP BY\n" + " `.theta`, `.sparse_sim`, `.dense_sim`\n" + ")"; + + query5.parse(); + for (auto it=env.cbegin(); it!=env.cend(); it++) + query5.template_defaults[it->first.c_str()] = it->second.c_str(); + + query5.template_defaults["table"] = "OriginalMC"; + query5.template_defaults["column"] = "DEGREES(Theta)"; + query5.template_defaults["sparse"] = str_sparse.c_str(); + query5.template_defaults["dense"] = str_dense.c_str(); + query5.template_defaults["theta"] = str_theta.c_str(); + query5.template_defaults["spectrum"] = spectrum.c_str(); + + if (print_queries) + PrintQuery(query5.str()); + + qlog << query5 << ";\n" << endl; + if (connection.connected()) + { + cout << query5.execute().info() << endl; + ShowWarnings(connection); + Dump(flog, connection, "SummaryOriginalMC"); + + const auto sec5 = Time().UnixTime()-start5.UnixTime(); + cout << "Execution time: " << sec5 << "s\n\n"; + } + + // ------------------------------------------------------------------- + + cout << separator("SummaryEventsMC") << '\n'; + + Time start5b; + + query5.template_defaults["table"] = "EventsMC"; + query5.template_defaults["column"] = "Zd"; + + if (print_queries) + PrintQuery(query5.str()); + + qlog << query5 << ";\n" << endl; + if (connection.connected()) + { + cout << query5.execute().info() << endl; + ShowWarnings(connection); + Dump(flog, connection, "SummaryEventsMC"); + + const auto sec5b = Time().UnixTime()-start5b.UnixTime(); + cout << "Execution time: " << sec5b << "s\n\n"; + } + + // ------------------------------------------------------------------- + // ---------------------------- ThetDist ----------------------------- + // ------------------------------------------------------------------- + + Time start6; + + // This table combines the analysis results vs. Binning in Estimated Energy and Simulated Energy + mysqlpp::Query query6(&connection); + query6 << + "CREATE TEMPORARY TABLE ThetaDist\n" + "ENGINE=MEMORY\n" + "AS\n" + "(\n" + " WITH ThetaCount AS\n" + " (\n" + " SELECT\n" + " `.theta`,\n" + " SUM(CountN) AS CountN\n" + " FROM\n" + " SummaryOriginalMC\n" " GROUP BY\n" " `.theta`\n" " )\n" " SELECT\n" - " `.theta`, lo, hi,\n" + " `.theta`,\n" " CountN,\n" " SQRT(CountN) AS ErrCountN,\n" @@ -1034,5 +1095,5 @@ " ObservationTime\n" " LEFT JOIN\n" - " EventCount USING(`.theta`)\n" + " ThetaCount USING (`.theta`)\n" " LEFT JOIN\n" " BinningTheta ON `.theta`=bin\n" @@ -1041,39 +1102,20 @@ ")"; - query4.parse(); + query6.parse(); for (auto it=env.cbegin(); it!=env.cend(); it++) - query4.template_defaults[it->first.c_str()] = it->second.c_str(); - - query4.template_defaults["bins"] = strzd.c_str(); + query6.template_defaults[it->first.c_str()] = it->second.c_str(); if (print_queries) - PrintQuery(query4.str()); - - qlog << query4 << ";\n" << endl; + PrintQuery(query6.str()); + + qlog << query6 << ";\n" << endl; if (connection.connected()) { - cout << query4.execute().info() << endl; + cout << query6.execute().info() << endl; ShowWarnings(connection); - Dump(flog, connection, "ThetaDistribution"); - - const auto sec4 = Time().UnixTime()-start4.UnixTime(); - cout << "Execution time: " << sec4 << "s\n"; - - - if (verbose>0) - { - const mysqlpp::StoreQueryResult res4 = connection.query("SELECT * FROM ThetaDistribution").store(); - - cout << " Bin MC Counts OnTime ZdWeight\n"; - const auto bins = binning_theta.vec(); - for (auto ir=res4.cbegin(); ir!=res4.cend(); ir++) - { - const mysqlpp::Row &row = *ir; - - const uint32_t bin = row[".theta"]; - cout << setw(5) << bins[bin] << ": " << setw(10) << row["CountN"] << " " << setw(10) << row["OnTime"] << " " << setw(10) << row["ZdWeight"] << '\n'; - } - cout << endl; - } + Dump(flog, connection, "ThetaDist"); + + const auto sec6 = Time().UnixTime()-start6.UnixTime(); + cout << "Execution time: " << sec6 << "s\n\n"; } @@ -1084,5 +1126,5 @@ .binningx = binning_theta, .binningy = {}, - .table = "ThetaDistribution", + .table = "ThetaDist", .x = ".theta", .y = "", @@ -1101,5 +1143,5 @@ .binningx = binning_theta, .binningy = {}, - .table = "ThetaDistribution", + .table = "ThetaDist", .x = ".theta", .y = "", @@ -1118,5 +1160,5 @@ .binningx = binning_theta, .binningy = {}, - .table = "ThetaDistribution", + .table = "ThetaDist", .x = ".theta", .y = "", @@ -1129,27 +1171,555 @@ }); - - // ------------------------------------------------------------------- - // ------------------------- 5: AnalysisData ------------------------- - // ------------------------------------------------------------------- - - cout << separator("AnalysisData") << '\n'; - - Time start5; - - /* 02:analyze-data.sql */ - mysqlpp::Query query5(&connection); - sindent indent5(query5); - query5 << - "CREATE TEMPORARY TABLE AnalysisData\n" + // ------------------------------------------------------------------- + // --------------------------- WeightedMC ---------------------------- + // ------------------------------------------------------------------- + + Time start7; + + // This table combines the analysis results vs. Binning in Estimated Energy and Simulated Energy + mysqlpp::Query query7(&connection); + query7 << + "CREATE TEMPORARY TABLE Weighted%100:table (\n" + " `.theta` SMALLINT UNSIGNED NOT NULL,\n" + " `.sparse_sim` SMALLINT UNSIGNED NOT NULL,\n" + " `.dense_sim` SMALLINT UNSIGNED NOT NULL,\n" + " CountN DOUBLE NOT NULL,\n" + " SumW DOUBLE NOT NULL,\n" + " SumW2 DOUBLE NOT NULL,\n" + " INDEX (`.theta`) USING HASH,\n" + " INDEX (`.sparse_sim`) USING HASH,\n" + " INDEX (`.dense_sim`) USING HASH\n" + ")\n" + "ENGINE=MEMORY\n" + "AS\n" "(\n" - " `.energy` SMALLINT UNSIGNED NOT NULL,\n" - " `Signal` DOUBLE NOT NULL,\n" - " `Background` DOUBLE NOT NULL,\n" - " `Excess` DOUBLE NOT NULL,\n" - " `Significance` DOUBLE NOT NULL,\n" - " `ErrExcess` DOUBLE NOT NULL,\n" + " SELECT\n" + " `.theta`,\n" + " `.sparse_sim`,\n" + " `.dense_sim`,\n" + " S.CountN,\n" + " S.SumW*ZdWeight AS SumW,\n" + " S.SumW2*POW(ErrZdWeight, 2) AS SumW2\n" + " FROM\n" + " Summary%100:table S\n" + " INNER JOIN\n" + " ThetaDist USING (`.theta`)\n" + ")"; + + query7.parse(); + for (auto it=env.cbegin(); it!=env.cend(); it++) + query7.template_defaults[it->first.c_str()] = it->second.c_str(); + + query7.template_defaults["table"] = "OriginalMC"; + + if (print_queries) + PrintQuery(query7.str()); + + qlog << query7 << ";\n" << endl; + if (connection.connected()) + { + cout << query7.execute().info() << endl; + ShowWarnings(connection); + Dump(flog, connection, "WeightedOriginalMC"); + + const auto sec7 = Time().UnixTime()-start7.UnixTime(); + cout << "Execution time: " << sec7 << "s\n\n"; + } + + // ------------------------------------------------------------------- + + Time start7b; + + query7.template_defaults["table"] = "EventsMC"; + + if (print_queries) + PrintQuery(query7.str()); + + qlog << query7 << ";\n" << endl; + if (connection.connected()) + { + cout << query7.execute().info() << endl; + ShowWarnings(connection); + Dump(flog, connection, "WeightedEventsMC"); + + const auto sec7b = Time().UnixTime()-start7b.UnixTime(); + cout << "Execution time: " << sec7b << "s\n\n"; + } + + // ------------------------------------------------------------------- + // --------------------------- AnalysisMC ---------------------------- + // ------------------------------------------------------------------- + + cout << separator("AnalysisMC") << '\n'; + + Time start8; + + // This table combines the analysis results vs. Binning in Estimated Energy and Simulated Energy + mysqlpp::Query query8(&connection); + sindent indent8(query8); + query8 << +/* + "CREATE TEMPORARY TABLE AnalysisMC\n" + "(\n" + " `.sparse` SMALLINT UNSIGNED NOT NULL,\n" + " `.dense` SMALLINT UNSIGNED NOT NULL,\n" + " SignalW DOUBLE NOT NULL,\n" // vs Eest (for spectral analysis) + " SignalN DOUBLE NOT NULL,\n" // vs Eest + " BackgroundN DOUBLE NOT NULL,\n" // vs Eest + " BackgroundW DOUBLE NOT NULL,\n" // vs Eest + " ExcessW DOUBLE NOT NULL,\n" // vs Eest + " ExcessN DOUBLE NOT NULL,\n" // vs Eest + " ErrExcessW DOUBLE NOT NULL,\n" // vs Eest + " ErrExcessN DOUBLE NOT NULL,\n" // vs Eest + " ThresholdW DOUBLE NOT NULL,\n" // vs Esim (for simulation analysis) + " ThresholdN DOUBLE NOT NULL,\n" // vs Esim + " BiasEst DOUBLE NOT NULL,\n" // vs Eest + " BiasSim DOUBLE NOT NULL,\n" // vs Esim + " ResolutionEst DOUBLE,\n" + " ResolutionSim DOUBLE,\n" + " INDEX (`.sparse`) USING HASH,\n" + " INDEX (`.dense`) USING HASH\n" + ") ENGINE=Memory\n" +*/ + "CREATE TEMPORARY TABLE AnalysisMC ENGINE=MEMORY\n" + "AS\n" + "(\n" + " WITH Excess AS\n" + " (\n" << indent(6) + << ifstream(analysis_sql).rdbuf() << indent(0) << + " ),\n" << + " Result AS\n" + " (\n" << indent(6) + << simulation_sql << indent(0) << // Must end with EOL and not in the middle of a comment + " )\n" + " SELECT * FROM Result\n" + ")"; + + query8.parse(); + for (auto it=env.cbegin(); it!=env.cend(); it++) + query8.template_defaults[it->first.c_str()] = it->second.c_str(); + + //query6.template_defaults["columns"] = "FileId, EvtNumber, CorsikaNumReuse,"; + query8.template_defaults["columns"] = "Zd, Energy, SpectralIndex,"; + query8.template_defaults["files"] = "MonteCarloFiles"; + query8.template_defaults["runinfo"] = "factmc.RunInfoMC"; + query8.template_defaults["events"] = "factmc.EventsMC"; + query8.template_defaults["positions"] = "factmc.PositionMC"; + + query8.template_defaults["sparse"] = str_sparse.c_str(); + query8.template_defaults["dense"] = str_dense.c_str(); + query8.template_defaults["theta"] = str_theta.c_str(); + query8.template_defaults["spectrum"] = spectrum.c_str(); + query8.template_defaults["estimator"] = estimator.c_str(); + + if (print_queries) + PrintQuery(query8.str()); + + qlog << query8 << ";\n" << endl; + if (connection.connected()) + { + cout << query8.execute().info() << endl; + ShowWarnings(connection); + Dump(flog, connection, "AnalysisMC"); + + const auto sec8 = Time().UnixTime()-start8.UnixTime(); + cout << "Execution time: " << sec8 << "s\n\n"; + } + + + // ------------------------------------------------------------------- + // ------------------------- SimulatedSpectrum ----------------------- + // ------------------------------------------------------------------- + + const vector binnings = { "dense", "sparse", "theta" }; + + for (auto ib=binnings.cbegin(); ib!=binnings.cend(); ib++) + { + const string table = "Summary"+(*ib=="theta" ? "Theta" : "TrueEnergy_"+*ib); + + cout << separator(table) << '\n'; + + // [lo^(1+gammaa) - hi^(1+gamma)] / (1+gamma) + + Time start9; + + /* + "CREATE TEMPORARY TABLE SummarySimulated\n" + "(\n" + " `.energy` SMALLINT UNSIGNED NOT NULL COMMENT 'Bin Index [MC Energy]',\n" + " CountN DOUBLE NOT NULL,\n" // COMMENT 'Number of events',\n" + " CountW DOUBLE NOT NULL,\n" // COMMENT 'Sum of weights, reweighted sum',\n" + " CountW2 DOUBLE NOT NULL,\n" // COMMENT 'Sum of square of weights'\n" " PRIMARY KEY (`.energy`) USING HASH\n" ") ENGINE=Memory\n" + */ + + mysqlpp::Query query9(&connection); + sindent indent9(query9); + query9 << + "CREATE TEMPORARY TABLE %100:table ENGINE=MEMORY " + "AS\n" + "(\n" + << indent(3) << summary_sim_sql << indent(0) << + ")"; + + // [ Sa^2/a^2 + Sb^2/b^2 ] * a^2/b^2 + // [ (sc^2)/c^2+(sb^2)/b^2+(sa^2)/a^2 ] * a^2*b^2/c^2 + + // [ 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 + + query9.parse(); + for (auto it=env.cbegin(); it!=env.cend(); it++) + query9.template_defaults[it->first.c_str()] = it->second.c_str(); + + query9.template_defaults["table"] = table.c_str(); + query9.template_defaults["binning"] = *ib=="theta" ? "BinningTheta" : ("BinningEnergy_"+*ib).c_str(); + query9.template_defaults["bin"] = *ib=="theta" ? "`.theta`" : ("`."+*ib+"_sim`").c_str(); + query9.template_defaults["binwidth"] = *ib=="theta" ? "1" : "1000/(POW(10,hi)-POW(10,lo))"; + + if (print_queries) + PrintQuery(query9.str()); + + qlog << query9 << ";\n" << endl; + if (connection.connected()) + { + cout << query9.execute().info() << endl; + ShowWarnings(connection); + Dump(flog, connection, table); + + const auto sec9 = Time().UnixTime()-start9.UnixTime(); + cout << "Execution time: " << sec9 << "s\n"; + } + + Histogram hist_sim; + hist_sim.table = table; + hist_sim.dir = *ib=="theta" ? "MC/theta" : "MC/"+*ib+"/TrueEnergy"; + hist_sim.x = *ib=="theta" ? ".theta" : "."+*ib+"_sim"; + hist_sim.axisx = *ib=="theta" ? "Zenith Angle #theta [#circ]" : "Energy lg(E_{sim}/GeV)"; + hist_sim.stats = false; + + if (*ib=="theta") + hist_sim.binningx=binning_theta; + if (*ib=="dense") + hist_sim.binningx=binning_dense; + if (*ib=="sparse") + hist_sim.binningx=binning_sparse; + + hist_sim.axisy = "Counts"; + + hist_sim.title = ""; + hist_sim.v = "SimCountN"; + hist_sim.err = "ErrSimCountN"; + WriteHistogram(connection, hist_sim); + + hist_sim.title = ""; + hist_sim.v = "TrigCountN"; + hist_sim.err = "ErrTrigCountN"; + WriteHistogram(connection, hist_sim); + + hist_sim.title = ""; + hist_sim.v = "SignalN"; + hist_sim.err = "ErrSignalN"; + WriteHistogram(connection, hist_sim); + + + hist_sim.axisy = *ib=="theta"?"dN/dE [cm^{-2} s^{-1}]":"dN/dE [cm^{-2} s^{-1} TeV^{-1}]"; + + hist_sim.title = ""; + hist_sim.v = "SimFluxW"; + hist_sim.err = "ErrSimFluxW"; + WriteHistogram(connection, hist_sim); + + hist_sim.title = ""; + hist_sim.v = "TrigFluxW"; + hist_sim.err = "ErrTrigFluxW"; + WriteHistogram(connection, hist_sim); + + hist_sim.title = ""; + hist_sim.v = "ExcessFluxW"; + hist_sim.err = "ErrExcessFluxW"; + WriteHistogram(connection, hist_sim); + + hist_sim.title = ""; + hist_sim.v = "SignalFluxW"; + hist_sim.err = "ErrSignalFluxW"; + WriteHistogram(connection, hist_sim); + + hist_sim.title = ""; + hist_sim.v = "BackgroundFluxW"; + hist_sim.err = "ErrBackgroundFluxW"; + WriteHistogram(connection, hist_sim); + + + hist_sim.title = ""; + hist_sim.v = "AvgEnergySimW"; + hist_sim.err = ""; + hist_sim.axisy = "/GeV"; + WriteHistogram(connection, hist_sim); + + hist_sim.title = ""; + hist_sim.v = "AvgEnergyEstW"; + hist_sim.err = ""; + hist_sim.axisy = "/GeV"; + WriteHistogram(connection, hist_sim); + + + hist_sim.title = ""; + hist_sim.v = "EffectiveAreaN"; + hist_sim.err = "ErrEffectiveAreaN"; + hist_sim.axisy = "A_{eff} [cm^{2}]"; + WriteHistogram(connection, hist_sim); + + hist_sim.title = ""; + hist_sim.v = "EffectiveAreaW"; + hist_sim.err = "ErrEffectiveAreaW"; + hist_sim.axisy = "A_{eff} [cm^{2}]"; + WriteHistogram(connection, hist_sim); + + + hist_sim.title = ""; + hist_sim.v = "BiasW"; + hist_sim.err = "ErrBiasW"; + hist_sim.axisy = ""; + WriteHistogram(connection, hist_sim); + + hist_sim.title = ""; + hist_sim.v = "ResolutionW"; + hist_sim.err = ""; + hist_sim.axisy = "#sigma(lg(E_{est}/E_{sim}))"; + WriteHistogram(connection, hist_sim); + + hist_sim.title = ""; + hist_sim.v = "TriggerEfficiencyN"; + hist_sim.err = "ErrTriggerEfficiencyN"; + hist_sim.axisy = "Efficiency"; + WriteHistogram(connection, hist_sim); + + hist_sim.title = ""; + hist_sim.v = "TriggerEfficiencyW"; + hist_sim.err = "ErrTriggerEfficiencyW"; + hist_sim.axisy = "Efficiency"; + WriteHistogram(connection, hist_sim); + + hist_sim.title = ""; + hist_sim.v = "CutEfficiencyN"; + hist_sim.err = "ErrCutEfficiencyN"; + hist_sim.axisy = "Efficiency"; + WriteHistogram(connection, hist_sim); + + hist_sim.title = ""; + hist_sim.v = "CutEfficiencyW"; + hist_sim.err = "ErrCutEfficiencyW"; + hist_sim.axisy = "Efficiency"; + WriteHistogram(connection, hist_sim); + + if (*ib=="theta") + continue; + + // ------------------------------------------------------------------- + // ------------------------- SimulatedSpectrum ----------------------- + // ------------------------------------------------------------------- + + cout << separator("SummaryEstimatedEnergy_"+*ib) << '\n'; + + // [lo^(1+gammaa) - hi^(1+gamma)] / (1+gamma) + + Time start10; + + mysqlpp::Query query10(&connection); + sindent indent10(query10); + query10 << + "CREATE TEMPORARY TABLE SummaryEstimatedEnergy_%100:binning ENGINE=MEMORY " + "AS\n" + "(\n" + << indent(3) << summary_est_sql << indent(6) << + ")"; + + query10.parse(); + for (auto it=env.cbegin(); it!=env.cend(); it++) + query10.template_defaults[it->first.c_str()] = it->second.c_str(); + + query10.template_defaults["binning"] = ib->c_str(); + + if (print_queries) + PrintQuery(query10.str()); + + qlog << query10 << ";\n" << endl; + if (connection.connected()) + { + cout << query10.execute().info() << endl; + ShowWarnings(connection); + Dump(flog, connection, "SummaryEstimatedEnergy_"+*ib); + + const auto sec10 = Time().UnixTime()-start10.UnixTime(); + cout << "Execution time: " << sec10 << "s\n"; + } + + Histogram hist_est; + hist_est.dir = "MC/"+*ib+"/EstimatedEnergy"; + hist_est.binningx = *ib=="dense"?binning_dense:binning_sparse; + hist_est.table = "SummaryEstimatedEnergy_"+*ib; + hist_est.x = "."+*ib+"_est"; + hist_est.axisx = "Energy lg(E_{est}/GeV)"; + hist_est.stats = false; + + hist_est.axisy = "Counts"; + + hist_est.title = ""; + hist_est.v = "SignalN"; + hist_est.err = "ErrSignalN"; + WriteHistogram(connection, hist_est); + + hist_est.title = ""; + hist_est.v = "BackgroundN"; + hist_est.err = "ErrBackgroundN"; + WriteHistogram(connection, hist_est); + + hist_est.title = ""; + hist_est.v = "ExcessN"; + hist_est.err = "ErrExcessN"; + WriteHistogram(connection, hist_est); + + + hist_est.title = ""; + hist_est.v = "AvgEnergySimW"; + hist_est.err = ""; + hist_est.axisy = "/GeV"; + WriteHistogram(connection, hist_est); + + hist_est.title = ""; + hist_est.v = "AvgEnergyEstW"; + hist_est.err = ""; + hist_est.axisy = "/GeV"; + WriteHistogram(connection, hist_est); + + + hist_est.axisy = "dN/dE [cm^{-2} s^{-1} TeV^{-1}]"; + + hist_est.title = ""; + hist_est.v = "SignalFluxW"; + hist_est.err = "ErrSignalFluxW"; + WriteHistogram(connection, hist_est); + + hist_est.title = ""; + hist_est.v = "BackgroundFluxW"; + hist_est.err = "ErrBackgroundFluxW"; + WriteHistogram(connection, hist_est); + + hist_est.title = ""; + hist_est.v = "ExcessFluxW"; + hist_est.err = "ErrExcessFluxW"; + WriteHistogram(connection, hist_est); + + + hist_est.title = ""; + hist_est.v = "BiasW"; + hist_est.err = "ErrBiasW"; + hist_est.axisy = ""; + WriteHistogram(connection, hist_est); + + hist_est.title = ""; + hist_est.v = "ResolutionW"; + hist_est.err = ""; + hist_est.axisy = "#sigma(lg(E_{est}/E_{sim}))"; + WriteHistogram(connection, hist_est); + + + // ------------------------------------------------------------------- + // ------------------------- SimulatedSpectrum ----------------------- + // ------------------------------------------------------------------- + + cout << separator("EnergyMigration_"+*ib) << '\n'; + + Time start11; + + mysqlpp::Query query11(&connection); + query11 << + "CREATE TEMPORARY TABLE EnergyMigration_%100:binning ENGINE=MEMORY\n" + "AS\n" + "(\n" + " SELECT\n" + " `.%100:binning:_est`,\n" + " `.%100:binning:_sim`,\n" + " SUM(SignalN) AS SignalN\n" + " FROM\n" + " AnalysisMC\n" + " GROUP BY\n" + " `.%100:binning:_est`, `.%100:binning:_sim`\n" + " ORDER BY\n" + " `.%100:binning:_est`, `.%100:binning:_sim`\n" + ")"; + + query11.parse(); + for (auto it=env.cbegin(); it!=env.cend(); it++) + query11.template_defaults[it->first.c_str()] = it->second.c_str(); + + query11.template_defaults["binning"] = ib->c_str(); + + if (print_queries) + PrintQuery(query11.str()); + + qlog << query11 << ";\n" << endl; + if (connection.connected()) + { + cout << query11.execute().info() << endl; + ShowWarnings(connection); + Dump(flog, connection, "EnergyMigration_"+*ib); + + const auto sec11 = Time().UnixTime()-start11.UnixTime(); + cout << "Execution time: " << sec11 << "s\n"; + } + + WriteHistogram(connection, { + .name = "Migration", + .title = "", + .dir = "MC/"+*ib, + .binningx = *ib=="dense"?binning_dense:binning_sparse, + .binningy = *ib=="dense"?binning_dense:binning_sparse, + .table = "EnergyMigration_"+*ib, + .x = "."+*ib+"_sim", + .y = "."+*ib+"_est", + .v = "SignalN", + .err = "", + .axisx = "Energy lg(E_{sim}/GeV)", + .axisy = "Energy lg(E_{est}/GeV)", + .axisz = "Counts", + .stats = false + }); + } + + if (mc_only) + { + cout << separator("Summary") << '\n'; + const auto sec = Time().UnixTime()-start.UnixTime(); + cout << "Total execution time: " << sec << "s\n" << endl; + + return 0; + } + + // ------------------------------------------------------------------- + // --------------------------- SummaryData --------------------------- + // ------------------------------------------------------------------- + + cout << separator("SummaryData") << '\n'; + + Time start12; + + mysqlpp::Query query12(&connection); + sindent indent12(query12); + query12 << + "CREATE TEMPORARY TABLE SummaryData ENGINE=MEMORY\n" +/* "(\n" + " `.theta` SMALLINT UNSIGNED NOT NULL,\n" + " `.sparse_est` SMALLINT UNSIGNED NOT NULL,\n" + " `Signal` DOUBLE NOT NULL,\n" + " `ErrSignal` DOUBLE NOT NULL,\n" + " `Background` DOUBLE NOT NULL,\n" + " `ErrBackground` DOUBLE NOT NULL,\n" + " `Excess` DOUBLE NOT NULL,\n" + " `ErrExcess` DOUBLE NOT NULL,\n" + " `Significance` DOUBLE NOT NULL,\n" + " PRIMARY KEY (`.sparse_est`) USING HASH\n" + ") ENGINE=Memory\n"*/ "AS\n" "(\n" @@ -1160,298 +1730,54 @@ " Result AS\n" " (\n" << indent(6) - << ifstream(data_sql).rdbuf() << indent(0) << // Must end with EOL and not in the middle of a comment + << data_sql << indent(0) << // Must end with EOL and not in the middle of a comment " )\n" " SELECT * FROM Result\n" ")"; - query5.parse(); + query12.parse(); for (auto it=env.cbegin(); it!=env.cend(); it++) - query5.template_defaults[it->first.c_str()] = it->second.c_str(); + query12.template_defaults[it->first.c_str()] = it->second.c_str(); //query5.template_defaults["columns"] = "FileId, EvtNumber,"; - query5.template_defaults["columns"] = ""; - query5.template_defaults["files"] = "DataFiles"; - query5.template_defaults["runinfo"] = "factdata.RunInfo"; - query5.template_defaults["events"] = "factdata.Images"; - query5.template_defaults["positions"] = "factdata.Position"; - - query5.template_defaults["bins"] = stre.c_str(); - query5.template_defaults["estimator"] = estimator.c_str(); + query12.template_defaults["columns"] = "fZenithDistanceMean,"; + query12.template_defaults["files"] = "DataFiles"; + query12.template_defaults["runinfo"] = "factdata.RunInfo"; + query12.template_defaults["events"] = "factdata.Images"; + query12.template_defaults["positions"] = "factdata.Position"; + + query12.template_defaults["sparse"] = str_sparse.c_str(); + query12.template_defaults["theta"] = str_theta.c_str(); + query12.template_defaults["estimator"] = estimator.c_str(); if (print_queries) - PrintQuery(query5.str()); - - qlog << query5 << ";\n" << endl; + PrintQuery(query12.str()); + + qlog << query12 << ";\n" << endl; if (connection.connected()) { - cout << query5.execute().info() << endl; + cout << query12.execute().info() << endl; ShowWarnings(connection); - Dump(flog, connection, "AnalysisData"); - - const auto sec5 = Time().UnixTime()-start5.UnixTime(); - cout << "Execution time: " << sec5 << "s\n"; - - if (verbose>0) - { - const mysqlpp::StoreQueryResult res5 = connection.query("SELECT * FROM AnalysisData").store(); - - cout << " Bin Signal Background Excess Significance Error" << endl; - for (auto row=res5.cbegin(); row!=res5.cend(); row++) - { - for (auto it=row->begin(); it!=row->end(); it++) - cout << setw(10) << *it << " "; - cout << '\n'; - } - cout << endl; - } - } - - - // ------------------------------------------------------------------- - // ----------------------------- 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'; - - Time start6; - - /* 02:analyze-simulation.sql */ - - // This table combines the analysis results vs. Binning in Estimated Energy and Simulated Energy - mysqlpp::Query query6(&connection); - sindent indent6(query6); - query6 << - "CREATE TEMPORARY TABLE AnalysisMC\n" - "(\n" - " `.energyest` SMALLINT UNSIGNED NOT NULL,\n" - " `.energysim` SMALLINT UNSIGNED NOT NULL,\n" - " SignalW DOUBLE NOT NULL,\n" - " BackgroundW DOUBLE NOT NULL,\n" - " ThresholdW DOUBLE NOT NULL,\n" - " SignalN DOUBLE NOT NULL,\n" - " BackgroundN DOUBLE NOT NULL,\n" - " ThresholdN DOUBLE NOT NULL,\n" - " ExcessW DOUBLE NOT NULL,\n" - " ExcessN DOUBLE NOT NULL,\n" - " ErrExcessW DOUBLE NOT NULL,\n" - " ErrExcessN DOUBLE NOT NULL,\n" - " BiasEst DOUBLE NOT NULL,\n" - " BiasSim DOUBLE NOT NULL,\n" - " ResolutionEst DOUBLE,\n" - " ResolutionSim DOUBLE,\n" - " INDEX (`.energyest`) USING HASH,\n" - " INDEX (`.energysim`) USING HASH\n" - ") ENGINE=Memory\n" - "AS\n" - "(\n" - " WITH Excess AS\n" - " (\n" << indent(6) - << ifstream(analysis_sql).rdbuf() << indent(0) << - " ),\n" << - " Result AS\n" - " (\n" << indent(6) - << ifstream(simulation_sql).rdbuf() << indent(0) << // Must end with EOL and not in the middle of a comment - " )\n" - " SELECT * FROM Result\n" - ")"; - - query6.parse(); - for (auto it=env.cbegin(); it!=env.cend(); it++) - query6.template_defaults[it->first.c_str()] = it->second.c_str(); - - //query6.template_defaults["columns"] = "FileId, EvtNumber, CorsikaNumReuse,"; - query6.template_defaults["columns"] = "Zd, Energy, SpectralIndex,"; - query6.template_defaults["files"] = "MonteCarloFiles"; - query6.template_defaults["runinfo"] = "factmc.RunInfoMC"; - query6.template_defaults["events"] = "factmc.EventsMC"; - query6.template_defaults["positions"] = "factmc.PositionMC"; - - query6.template_defaults["energyest"] = stre.c_str(); - query6.template_defaults["energysim"] = strth.c_str(); - query6.template_defaults["theta"] = strzd.c_str(); - query6.template_defaults["spectrum"] = spectrum.c_str(); - query6.template_defaults["estimator"] = estimator.c_str(); - - if (print_queries) - PrintQuery(query6.str()); - - qlog << query6 << ";\n" << endl; - if (connection.connected()) - { - cout << query6.execute().info() << endl; - ShowWarnings(connection); - Dump(flog, connection, "AnalysisMC"); - - const auto sec6 = Time().UnixTime()-start6.UnixTime(); - cout << "Execution time: " << sec6 << "s\n\n"; - } - - // ------------------------------------------------------------------- - // ------------------------- SimulatedSpectrum ----------------------- - // ------------------------------------------------------------------- - - cout << separator("SimulatedSpectrum") << '\n'; - - // [lo^(1+gammaa) - hi^(1+gamma)] / (1+gamma) - - Time start7; - /* 02:get-corsika-events.sql */ - - // FIXME: Theta weights? - // FIXME: energysim binning - mysqlpp::Query query7(&connection); - query7 << - "CREATE TEMPORARY TABLE SimulatedSpectrum\n" - "(\n" - " `.energy` SMALLINT UNSIGNED NOT NULL COMMENT 'Bin Index [MC Energy]',\n" - " CountN DOUBLE NOT NULL,\n" // COMMENT 'Number of events',\n" - " CountW DOUBLE NOT NULL,\n" // COMMENT 'Sum of weights, reweighted sum',\n" - " CountW2 DOUBLE NOT NULL,\n" // COMMENT 'Sum of square of weights'\n" - " PRIMARY KEY (`.energy`) USING HASH\n" - ") ENGINE=Memory\n" - "AS\n" - "(\n" - " SELECT\n" - " INTERVAL(LOG10(Energy), %100:energyest) AS `.energy`,\n" - " COUNT(*) AS CountN,\n" - " SUM( (%101:spectrum)/pow(Energy, SpectralIndex) ) AS CountW,\n" // Contents is: CountW - " SUM(POW((%101:spectrum)/pow(Energy, SpectralIndex),2)) AS CountW2\n" // Error is: SQRT(CountW2) - " FROM\n" - " MonteCarloFiles\n" - " LEFT JOIN\n" - " factmc.RunInfoMC USING (FileId)\n" - " LEFT JOIN\n" - " factmc.OriginalMC USING (FileId)\n" - " GROUP BY\n" - " `.energy`\n" - " ORDER BY\n" - " `.energy`\n" - ")"; - - query7.parse(); - for (auto it=env.cbegin(); it!=env.cend(); it++) - query7.template_defaults[it->first.c_str()] = it->second.c_str(); - - //query7.template_defaults["list"] = listmc.c_str(); - query7.template_defaults["energyest"] = stre.c_str(); - query7.template_defaults["spectrum"] = spectrum.c_str(); - - if (print_queries) - PrintQuery(query7.str()); - - qlog << query7 << ";\n" << endl; - if (connection.connected()) - { - cout << query7.execute().info() << endl; - ShowWarnings(connection); - Dump(flog, connection, "SimulatedSpectrum"); - - const auto sec7 = Time().UnixTime()-start7.UnixTime(); - cout << "Execution time: " << sec7 << "s\n"; - - - if (verbose>0) - { - const mysqlpp::StoreQueryResult res7 = connection.query("SELECT * FROM SimulatedSpectrum").store(); - - cout << " Bin CountW CountW2" << endl; - const auto bins = binning_esim.vec(); - for (auto ir=res7.cbegin(); ir!=res7.cend(); ir++) - { - const mysqlpp::Row &row = *ir; - - const uint32_t bin = row[".energy"]; - cout << setw(5) << bins[bin] << ": " << setw(10) << row["CountW"] << " " << setw(10) << row["CountW2"] << '\n'; - } - cout << endl; - } - } - - - // ------------------------------------------------------------------- - // ----------------------------- Spectrum ---------------------------- - // ------------------------------------------------------------------- - - cout << separator("Spectrum") << '\n'; - - Time start8; - - /* 02:calculate-spectrum.sql */ - - mysqlpp::Query query8(&connection); - query8 << - // SELECT SUM(CountN) AS CountN, SUM(OnTime) AS OnTime FROM ThetaHist - "CREATE TEMPORARY TABLE Spectrum\n" + Dump(flog, connection, "SummaryData"); + + const auto sec12 = Time().UnixTime()-start12.UnixTime(); + cout << "Execution time: " << sec12 << "s\n"; + } + + // ------------------------------------------------------------------- + // --------------------------- Spectrum ------------------------------ + // ------------------------------------------------------------------- + + const vector targets = { "Theta", "Energy" }; + + for (auto ib=targets.cbegin(); ib!=targets.cend(); ib++) + { + const string table = "Spectrum"+*ib; + + cout << separator(table) << '\n'; + + Time start13; + + /* + "CREATE TEMPORARY TABLE Spectrum\n" "(\n" " `.energy` SMALLINT UNSIGNED NOT NULL COMMENT 'Bin Index [Energy]' PRIMARY KEY,\n" @@ -1473,741 +1799,201 @@ " ErrSignalW DOUBLE NOT NULL COMMENT 'Error of weighted number of signal events in simulated data',\n" " ErrBackgroundW DOUBLE NOT NULL COMMENT 'Error of weighted number of background events in simulated data',\n" - " Flux DOUBLE NOT NULL COMMENT 'dN/dA/dt [cm^-2 s-^1 TeV^-1]',\n" - " ErrFlux DOUBLE NOT NULL COMMENT 'dN/dA/dt [cm^-2 s-^1 TeV^-1]',\n" + " Flux DOUBLE NOT NULL COMMENT 'Measured Monte Carlo Flux dN/dA/dt [cm^-2 s-^1 TeV^-1]',\n" + " ErrFlux DOUBLE NOT NULL COMMENT 'Error of measures Monte Carlo Flux dN/dA/dt [cm^-2 s-^1 TeV^-1]',\n" + //" CountSim DOUBLE NOT NULL COMMENT 'Simulated Monte Carlo Events',\n" + //" ErrCountSim DOUBLE NOT NULL COMMENT 'Error of Simulated Monte Carlo Events',\n" + //" FluxSim DOUBLE NOT NULL COMMENT 'Simulated Monte Carlo Flux dN/dA/dt [cm^-2 s-^1 TeV^-1]',\n" + //" ErrFluxSim DOUBLE NOT NULL COMMENT 'Error of Simulated Monte Carlo Flux dN/dA/dt [cm^-2 s-^1 TeV^-1]',\n" " Bias DOUBLE NOT NULL COMMENT 'Energy Bias, average residual in lg(E)',\n" " Resolution DOUBLE NOT NULL COMMENT 'Energy resolution, standard divation of residual in lg(E)',\n" - " EfficiencyN DOUBLE NOT NULL COMMENT 'Simulated cut efficiency (weighted)',\n" - " EfficiencyW DOUBLE NOT NULL COMMENT 'Simulated cut efficiency (unweighted)',\n" - " ErrEfficiencyN DOUBLE NOT NULL COMMENT 'Error of simulated cut efficiency (weighted)',\n" - " ErrEfficiencyW DOUBLE NOT NULL COMMENT 'Error of simulated cut efficiency (unweighted)'\n" + //" EfficiencyN DOUBLE NOT NULL COMMENT 'Simulated cut efficiency (weighted)',\n" + //" EfficiencyW DOUBLE NOT NULL COMMENT 'Simulated cut efficiency (unweighted)',\n" + //" ErrEfficiencyN DOUBLE NOT NULL COMMENT 'Error of simulated cut efficiency (weighted)',\n" + //" ErrEfficiencyW DOUBLE NOT NULL COMMENT 'Error of simulated cut efficiency (unweighted)'\n" ") ENGINE=Memory\n" - "AS\n" - "(\n" - " WITH ThetaSums AS\n" - " (\n" - " SELECT\n" - " SUM(CountN) AS CountSim,\n" - " SUM(OnTime) AS ObsTime\n" - " FROM\n" - " ThetaDistribution\n" - " ),\n" - " Area AS\n" - " (\n" - " SELECT\n" - " POW(MinImpactHi,2)*PI() AS Area\n" - " FROM\n" - " MonteCarloArea\n" - " ),\n" - " ResultMC AS\n" // Group AnalysisMC by EnergyEst Binning - " (\n" - " SELECT\n" - " `.energyest` AS `.energy`,\n" - " ANY_VALUE(SignalW) AS SignalW,\n" - " ANY_VALUE(SignalW2) AS SignalW2,\n" - " ANY_VALUE(BackgroundW) AS BackgroundW,\n" - " ANY_VALUE(BackgroundW2) AS BackgroundW2,\n" - " ANY_VALUE(SignalN) AS SignalN,\n" - " ANY_VALUE(BackgroundN) AS BackgroundN,\n" - " ANY_VALUE(ExcessW) AS ExcessW,\n" - " ANY_VALUE(ExcessN) AS ExcessN,\n" - " ANY_VALUE(ErrExcessW) AS ErrExcessW,\n" - " ANY_VALUE(ErrExcessN) AS ErrExcessN,\n" - " ANY_VALUE(BiasEst) AS Bias,\n" - " ANY_VALUE(ResolutionEst) AS Resolution\n" - " FROM\n" - " AnalysisMC\n" - " GROUP BY\n" - " `.energy`\n" - " ORDER BY\n" - " `.energy`\n" - " )\n" - " SELECT\n" - " `.energy`, lo, hi,\n" // Scale for Theta-Weights - " `Signal`, `Background`/5 AS `Background`, `Excess`, `ErrExcess`, `Significance`,\n" - " SQRT(`Signal`) AS ErrSignal,\n" - " SQRT(`SignalW2`) AS ErrSignalW,\n" - " SQRT(`Background`)/5 AS ErrBackground,\n" - " SQRT(`BackgroundW2`)/5 AS ErrBackgroundW,\n" - " ExcessN, ExcessW, ErrExcessN, ErrExcessW, SignalW, BackgroundW,\n" - " AnalysisData.Excess/ResultMC.ExcessW*SimulatedSpectrum.CountW * 1000/(POW(10,hi)-POW(10,lo)) /Area/ObsTime / CountSim*ObsTime AS Flux,\n" - " AnalysisData.Excess/ResultMC.ExcessW*SimulatedSpectrum.CountW * 1000/(POW(10,hi)-POW(10,lo)) /Area/ObsTime / CountSim*ObsTime\n" - " * SQRT(\n" - " + POW(AnalysisData.ErrExcess / AnalysisData.Excess, 2)\n" - " + POW(ResultMC.ErrExcessW / ResultMC.ExcessW, 2)\n" - " + SimulatedSpectrum.CountW2 / POW(SimulatedSpectrum.CountW,2)\n" - " ) AS ErrFlux,\n" - " Bias,\n" - " Resolution,\n" - " ResultMC.ExcessW/SimulatedSpectrum.CountW * CountSim/ObsTime AS EfficiencyW,\n" - " ResultMC.ExcessN/SimulatedSpectrum.CountN AS EfficiencyN,\n" - " ( POW(ResultMC.ErrExcessW/ResultMC.ExcessW, 2) + POW(SQRT(SimulatedSpectrum.CountW2)/SimulatedSpectrum.CountW, 2) )\n" - " * POW(ResultMC.ExcessW/SimulatedSpectrum.CountW, 2) * CountSim/ObsTime AS ErrEfficiencyW,\n" - " ( POW(ResultMC.ErrExcessN, 2) + POW(ResultMC.ExcessN, 2)/SimulatedSpectrum.CountN)/POW(SimulatedSpectrum.CountN, 2) AS ErrEfficiencyN,\n" - " 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" - " INNER JOIN\n" - " ResultMC USING(`.energy`)\n" - " INNER JOIN\n" - " SimulatedSpectrum USING(`.energy`)\n" - " INNER JOIN\n" - " BinningEnergyEst ON `.energy`=bin\n" - " CROSS JOIN\n" - " ThetaSums, Area\n" - " WHERE\n" - " AnalysisData.Excess>0\n" - " ORDER BY\n" - " `.energy`\n" - ")" - ; - - // [ Sa^2/a^2 + Sb^2/b^2 ] * a^2/b^2 - // [ (sc^2)/c^2+(sb^2)/b^2+(sa^2)/a^2 ] * a^2*b^2/c^2 - - query8.parse(); - for (auto it=env.cbegin(); it!=env.cend(); it++) - query8.template_defaults[it->first.c_str()] = it->second.c_str(); - - //query8.template_defaults["area"] = area; - //query8.template_defaults["ontime"] = resX[0]["OnTime"].data(); - //query8.template_defaults["count"] = resX[0]["CountN"].data(); - - if (print_queries) - PrintQuery(query8.str()); - - qlog << query8 << ";\n" << endl; - if (connection.connected()) - { - cout << query8.execute().info() << endl; - ShowWarnings(connection); - Dump(flog, connection, "Spectrum"); - - const auto sec8 = Time().UnixTime()-start8.UnixTime(); - cout << "Execution time: " << sec8 << "s\n"; - - - const mysqlpp::StoreQueryResult res8 = connection.query("SELECT * FROM Spectrum").store(); - - if (verbose>0) +*/ + + mysqlpp::Query query13(&connection); + sindent indent13(query13); + query13 << + "CREATE TEMPORARY TABLE %100:table ENGINE=MEMORY " + "AS\n" + "(\n" + << indent(3) << spectrum_sql << indent(0) << + ")"; +*/ + // [ Sa^2/a^2 + Sb^2/b^2 ] * a^2/b^2 + // [ (sc^2)/c^2+(sb^2)/b^2+(sa^2)/a^2 ] * a^2*b^2/c^2 + + query13.parse(); + for (auto it=env.cbegin(); it!=env.cend(); it++) + query13.template_defaults[it->first.c_str()] = it->second.c_str(); + + query13.template_defaults["table"] = table.c_str(); + query13.template_defaults["bin"] = *ib=="Theta" ? "`.theta`" : "`.sparse_est`"; + query13.template_defaults["id"] = *ib=="Theta" ? "Sim" : "Est"; + query13.template_defaults["weight"] = *ib=="Theta" ? "ZdWeight" : "1"; + query13.template_defaults["errweight"] = *ib=="Theta" ? "ErrZdWeight" : "1"; + if (*ib=="Theta") { - cout << " Bin Flux Error" << endl; - const auto bins = binning_eest.vec(); - for (auto ir=res8.cbegin(); ir!=res8.cend(); ir++) + query13.template_defaults["join1"] = "SummaryTheta Sim USING (`.theta`)"; + query13.template_defaults["join2"] = "ThetaDist USING (`.theta`)"; + } + else + { + query13.template_defaults["join1"] = "SummaryEstimatedEnergy_sparse Est USING (`.sparse_est`)"; + query13.template_defaults["join2"] = "SummaryTrueEnergy_sparse Sim ON Est.`.sparse_est`=Sim.`.sparse_sim`"; + } + + if (print_queries) + PrintQuery(query13.str()); + + qlog << query13 << ";\n" << endl; + if (connection.connected()) + { + cout << query13.execute().info() << endl; + ShowWarnings(connection); + Dump(flog, connection, table); + + const auto sec13 = Time().UnixTime()-start13.UnixTime(); + cout << "Execution time: " << sec13 << "s\n"; + + + const mysqlpp::StoreQueryResult res13 = connection.query("SELECT * FROM "+table).store(); + + // -------------------------------------------------------------------------- +#ifdef HAVE_ROOT + TFeldmanCousins fc; + fc.SetCL(confidence); + fc.SetMuStep(0.2); + // f.SetMuMax(10*(sig+bg)); //has to be higher than sig+bg!! + // Double_t ul=f.CalculateUpperLimit(x,y); + // x=Dat.Signal : number of observed events in the experiment + // y=Dat.Background/5 : average number of observed events in background region + + TRolke rolke(confidence); + // rolke.SetPoissonBkgBinomEff(x,y,z,tau,m) + // x=Dat.Signal : number of observed events in the experiment + // y=Dat.Background : number of observed events in background region + // tau=0.2 : the ratio between signal and background region + // m=Sim.SimFluxN : number of MC events generated + // z=Sim.AnaFluxN : number of MC events observed +#endif + // -------------------------------------------------------------------------- + + // Crab Nebula: 1 TeV: 3e-7 / m^2 / s / TeV + // Crab Nebula: 1 TeV: 3e-11 / cm^2 / s / TeV + + map feldman_ul; + map rolke_ul; + map rolke_ll; + + for (auto ir=res13.cbegin(); ir!=res13.cend(); ir++) { + // This is not efficient but easier to read and safer const mysqlpp::Row &row = *ir; - const uint32_t bin = row[".energy"]; - cout << setw(5) << bins[bin] << ": " << setw(10) << row["Flux"] << " " << setw(10) << row["ErrFlux"] << '\n'; + const size_t bin = row[*ib=="Theta" ? ".theta" : ".sparse_est"]; + const double flux = row["Flux"]; + const double error = row["ErrFlux"]; + +#ifdef HAVE_ROOT + const double dat_sig = row["Signal"]; + const double dat_bg = row["Background"]; + + const double eff = row["Efficiency"]; + const double scale = row["Scale"]; + + fc.SetMuMax(10*(dat_sig+dat_bg/5)); //has to be higher than sig+bg!! + rolke.SetPoissonBkgKnownEff(dat_sig, dat_bg, 5, eff); + + if (feldman) + feldman_ul[bin] = fc.CalculateUpperLimit(dat_sig, dat_bg/5)*scale/eff; + + rolke_ll[bin] = rolke.GetLowerLimit()*scale; + rolke_ul[bin] = rolke.GetUpperLimit()*scale; +#endif + if (verbose>0) + { + cout << setw(5) << row["center"] << ": "; + cout << setw(10) << row["Excess"] << " "; + cout << setw(10) << flux << " "; + cout << setw(10) << error << " "; + cout << setw(10) << row["Significance"] << '\n'; + } } - cout << endl; + + Histogram hist_zd; + hist_zd.dir = "Data/"+*ib; + hist_zd.table = table; + hist_zd.binningx = *ib=="Theta" ? binning_theta : binning_sparse; + hist_zd.x = *ib=="Theta" ? ".theta" : ".sparse_est"; + hist_zd.axisx = *ib=="Theta" ? "Zenith Distance #theta [#circ]" : "Energy lg(E/GeV)"; + hist_zd.stats = false; + + hist_zd.title = ""; + hist_zd.v = "Signal"; + hist_zd.err = "ErrSignal"; + hist_zd.axisy = "Counts"; + WriteHistogram(connection, hist_zd); + + hist_zd.title = ""; + hist_zd.v = "Background"; + hist_zd.err = "ErrBackground"; + hist_zd.axisy = "Counts"; + WriteHistogram(connection, hist_zd); + + hist_zd.title = ""; + hist_zd.v = "Excess"; + hist_zd.err = "ErrExcess"; + hist_zd.axisy = "Counts"; + WriteHistogram(connection, hist_zd); + + hist_zd.title = ""; + hist_zd.v = "Significance"; + hist_zd.err = ""; + hist_zd.axisy = "#sigma"; + WriteHistogram(connection, hist_zd); + + hist_zd.title = ""; + hist_zd.v = "AvgEnergyEst"; + hist_zd.err = ""; + hist_zd.axisy = "/GeV"; + WriteHistogram(connection, hist_zd); + + hist_zd.v = "ExcessRatio"; + hist_zd.err = "ErrExcessRatio"; + hist_zd.axisy = "Ratio"; + WriteHistogram(connection, hist_zd); + + + hist_zd.axisy = *ib=="Theta" ? "dN/dE [cm^{-2} s^{-1}]" : "dN/dE [cm^{-2} s^{-1} TeV^{-1}]"; + + hist_zd.name = "Spectrum"; + hist_zd.v = "Flux"; + hist_zd.err = "ErrFlux"; + WriteHistogram(connection, hist_zd); + +#ifdef HAVE_ROOT + hist_zd.axisy = *ib=="Theta" ? "UL [cm^{-2} s^{-1}]" : "UL [cm^{-2} s^{-1} TeV^{-1}]"; + + if (feldman) + { + hist_zd.name = "FeldmanCousins"; + WriteHistogram(connection, hist_zd, feldman_ul); + } + + hist_zd.name = "RolkeUL"; + WriteHistogram(connection, hist_zd, rolke_ul); + + hist_zd.name = "RolkeLL"; + WriteHistogram(connection, hist_zd, rolke_ll); +#endif } - - // -------------------------------------------------------------------------- - - // Crab Nebula: 1 TeV: 3e-7 / m^2 / s / TeV - // Crab Nebula: 1 TeV: 3e-11 / cm^2 / s / TeV - - sindent indentm(mlog); - - mlog << "void spectrum()\n"; - mlog << "{\n" << indent(4); - mlog << "// Energy Spectrum (e.g. Crab: 3e-11 [cm^-2 s-^1 TeV^-1] @ 1TeV)\n"; - mlog << "TGraphErrors g;\n"; - mlog << "g.SetNameTitle(\"Spectrum\", \"Energy Spectrum\");\n"; - for (auto ir=res8.cbegin(); ir!=res8.cend(); ir++) - { - // This is not efficient but easier to read and safer - const mysqlpp::Row &row = *ir; - - const double hi = row["hi"]; - const double lo = row["lo"]; - const double center = (hi+lo)/2; - const double flux = row["Flux"]; - const double error = row["ErrFlux"]; - - mlog << "g.SetPoint(g.GetN(), " << setw(8) << center << ", " << flux << ");\n"; - mlog << "g.SetPointError(g.GetN()-1, 0, " << error << ");\n"; - } - mlog << "g.GetXaxis()->SetTitle(\"lg(E/TeV)\");\n"; - mlog << "g.GetYaxis()->SetTitle(\"dN/dE [cm^{-2} s^{-1} TeV^{-1}]\");"; - mlog << endl; - } - - WriteHistogram(connection, { - .name = "Signal", - .title = "Signal", - .dir = "Eest/Measurement", - .binningx = binning_eest, - .binningy = {}, - .table = "Spectrum", - .x = ".energy", - .y = "", - .v = "Signal", - .err = "ErrSignal", - .axisx = "Energy lg(E_{est}/GeV)", - .axisy = "Counts", - .axisz = "", - .stats = true - }); - - WriteHistogram(connection, { - .name = "Background", - .title = "Background", - .dir = "Eest/Measurement", - .binningx = binning_eest, - .binningy = {}, - .table = "Spectrum", - .x = ".energy", - .y = "", - .v = "Background", - .err = "ErrBackground", - .axisx = "Energy lg(E_{est}/GeV)", - .axisy = "Count", - .axisz = "", - .stats = true - }); - - WriteHistogram(connection, { - .name = "Excess", - .title = "Excess", - .dir = "Eest/Measurement", - .binningx = binning_eest, - .binningy = {}, - .table = "Spectrum", - .x = ".energy", - .y = "", - .v = "Excess", - .err = "ErrExcess", - .axisx = "Energy lg(E_{est}/GeV)", - .axisy = "Signal - Background (Counts)", - .axisz = "", - .stats = true - }); - - WriteHistogram(connection, { - .name = "SignalW", - .title = "SignalW", - .dir = "Eest/Simulation/Weighted", - .binningx = binning_eest, - .binningy = {}, - .table = "Spectrum", - .x = ".energy", - .y = "", - .v = "SignalW", - .err = "ErrSignalW", - .axisx = "Energy lg(E_{est}/GeV)", - .axisy = "Weighted", - .axisz = "", - .stats = true - }); - - WriteHistogram(connection, { - .name = "BackgroundW", - .title = "BackgroundW", - .dir = "Eest/Simulation/Weighted", - .binningx = binning_eest, - .binningy = {}, - .table = "Spectrum", - .x = ".energy", - .y = "", - .v = "BackgroundW", - .err = "ErrBackgroundW", - .axisx = "Energy lg(E_{est}/GeV)", - .axisy = "Weighted", - .axisz = "", - .stats = true - }); - - WriteHistogram(connection, { - .name = "ExcessW", - .title = "ExcessW", - .dir = "Eest/Simulation/Weighted", - .binningx = binning_eest, - .binningy = {}, - .table = "Spectrum", - .x = ".energy", - .y = "", - .v = "ExcessW", - .err = "ErrExcessW", - .axisx = "Energy lg(E_{est}/GeV)", - .axisy = "Signal - Background (Weighted)", - .axisz = "", - .stats = true - }); - - WriteHistogram(connection, { - .name = "Significance", - .title = "Significance", - .dir = "Eest/Measurement", - .binningx = binning_eest, - .binningy = {}, - .table = "Spectrum", - .x = ".energy", - .y = "", - .v = "Significance", - .err = "", - .axisx = "Energy lg(E_{est}/GeV)", - .axisy = "Li/Ma Significance", - .axisz = "", - .stats = true - }); - - WriteHistogram(connection, { - .name = "Bias", - .title = "Energy Bias", - .dir = "Eest", - .binningx = binning_eest, - .binningy = {}, - .table = "Spectrum", - .x = ".energy", - .y = "", - .v = "Bias", - .err = "Resolution", - .axisx = "Energy lg(E_{sim}/GeV)", - .axisy = "lg(E_{est}/E_{sim})", - .axisz = "", - .stats = false - }); - - WriteHistogram(connection, { - .name = "Resolution", - .title = "Energy Resolution", - .dir = "Eest", - .binningx = binning_eest, - .binningy = {}, - .table = "Spectrum", - .x = ".energy", - .y = "", - .v = "Resolution", - .err = "", - .axisx = "Energy lg(E_{sim}/GeV)", - .axisy = "#sigma(lg(E_{est}/E_{sim}))", - .axisz = "", - .stats = false - }); - - WriteHistogram(connection, { - .name = "EfficiencyN", - .title = "Efficiency (Counts)", - .dir = "Eest/Efficiency", - .binningx = binning_eest, - .binningy = {}, - .table = "Spectrum", - .x = ".energy", - .y = "", - .v = "EfficiencyN", - .err = "ErrEfficiencyN", - .axisx = "Energy lg(E_{sim}/GeV)", - .axisy = "Efficiency", - .axisz = "", - .stats = true - }); - - WriteHistogram(connection, { - .name = "EfficiencyW", - .title = "Efficiency (Weighted)", - .dir = "Eest/Efficiency", - .binningx = binning_eest, - .binningy = {}, - .table = "Spectrum", - .x = ".energy", - .y = "", - .v = "EfficiencyW", - .err = "ErrEfficiencyW", - .axisx = "Energy lg(E_{sim}/GeV)", - .axisy = "Efficiency", - .axisz = "", - .stats = true - }); - - 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", - .dir = "", - .binningx = binning_eest, - .binningy = {}, - .table = "Spectrum", - .x = ".energy", - .y = "", - .v = "Flux", - .err = "ErrFlux", - .axisx = "Energy lg(E/GeV)", - .axisy = "dN/dE [cm^{-2} s^{-1} TeV^{-1}]", - .axisz = "", - .stats = false - }); - - - // ------------------------------------------------------------------- - // ---------------------------- Threshold ---------------------------- - // ------------------------------------------------------------------- - - cout << separator("Threshold") << '\n'; - - Time start9; - - /* 02:calculate-threshold.sql */ - - // This query gets the analysis results versus Simulated Energy from the combined table - mysqlpp::Query query9(&connection); - query9 << - // SELECT SUM(CountN) AS CountN, SUM(OnTime) AS OnTime FROM ThetaHist - "CREATE TEMPORARY TABLE Threshold ENGINE=Memory AS\n" - "(\n" - " WITH ThetaSums AS\n" - " (\n" - " SELECT\n" - " SUM(CountN) AS CountSim,\n" - " SUM(OnTime) AS ObsTime\n" - " FROM\n" - " ThetaDistribution\n" - " ),\n" - " Area AS\n" - " (\n" - " SELECT\n" - " POW(MinImpactHi,2)*PI() AS Area\n" - " FROM\n" - " MonteCarloArea\n" - " ),\n" - " ResultMC AS\n" // Group AnalysisMC by EnergySim Binning - " (\n" - " SELECT\n" - " `.energysim` AS `.energy`,\n" - " ANY_VALUE(ThresholdW) AS ThresholdW,\n" - " ANY_VALUE(ThresholdW2) AS ThresholdW2,\n" - " ANY_VALUE(ThresholdN) AS ThresholdN,\n" - " ANY_VALUE(BiasSim) AS Bias,\n" - " ANY_VALUE(ResolutionSim) AS Resolution\n" - " FROM\n" - " AnalysisMC\n" - " GROUP BY\n" - " `.energy`\n" - " )\n" - " SELECT\n" - " `.energy`, lo, hi,\n" - " ThresholdW,\n" - " SQRT(ThresholdW2) AS ErrThresholdW,\n" - " ThresholdN,\n" - " SQRT(ThresholdN) AS ErrThresholdN,\n" // Scale for Theta-Weights - " ThresholdW * 1000/(POW(10,hi)-POW(10,lo)) / Area / CountSim*ObsTime AS Flux,\n" - " SQRT(ThresholdW2) * 1000/(POW(10,hi)-POW(10,lo)) / Area / CountSim*ObsTime AS ErrFlux,\n" - " 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" - " CROSS JOIN\n" - " ThetaSums, Area\n" - " WHERE\n" - " ThresholdW>0 AND ThresholdW2>0\n" - " ORDER BY\n" - " `.energy`\n" - ")"; - - query9.parse(); - for (auto it=env.cbegin(); it!=env.cend(); it++) - query9.template_defaults[it->first.c_str()] = it->second.c_str(); - - //query9.template_defaults["area"] = area; - //query9.template_defaults["ontime"] = resX[0]["OnTime"].data(); - //query9.template_defaults["count"] = resX[0]["CountN"].data(); - - if (print_queries) - PrintQuery(query9.str()); - - qlog << query9 << ";\n" << endl; - if (connection.connected()) - { - cout << query9.execute().info() << endl; - ShowWarnings(connection); - Dump(flog, connection, "Threshold"); - - const auto sec9 = Time().UnixTime()-start9.UnixTime(); - cout << "Execution time: " << sec9 << "s\n"; - - // -------------------------------------------------------------------------- - - const mysqlpp::StoreQueryResult res9 = connection.query("SELECT * FROM Threshold").store(); - - sindent indentm(mlog, 4); - - mlog << '\n'; - mlog << "// Energy Threshold\n"; - mlog << "TGraphErrors g;\n"; - mlog << "g.SetNameTitle(\"Threshold\", \"Energy Threshold\");\n"; - for (auto ir=res9.cbegin(); ir!=res9.cend(); ir++) - { - // This is not efficient but easier to read and safer - const mysqlpp::Row &row = *ir; - - const double hi = row["hi"]; - const double lo = row["lo"]; - const double center = (hi+lo)/2; - const double width = pow(10, hi)-pow(10, lo); - const double flux = row["Flux"] /width; - const double error = row["ErrFlux"]/width; - - mlog << "g.SetPoint(g.GetN(), " << setw(8) << center << ", " << flux << ");\n"; - mlog << "g.SetPointError(g.GetN()-1, 0, " << error << ");\n"; - } - mlog << "g.GetXaxis()->SetTitle(\"lg(E/TeV)\");\n"; - mlog << "g.GetYaxis()->SetTitle(\"dN/dE [cm^{-2} s^{-1} TeV^{-1}]\");\n"; - mlog << indent(0) << "}" << endl; - } - - WriteHistogram(connection, { - .name = "SimExcessW", - .title = "Weighted Simulated Excess", - .dir = "Esim/Simulation/Weighted", - .binningx = binning_esim, - .binningy = {}, - .table = "Threshold", - .x = ".energy", - .y = "", - .v = "ThresholdW", - .err = "ErrThresholdW", - .axisx = "Energy lg(E_{sim}/GeV)", - .axisy = "Weighted Counts", - .axisz = "", - .stats = true - }); - - WriteHistogram(connection, { - .name = "SimExcessN", - .title = "Simulated Excess", - .dir = "Esim/Simulation/Counts", - .binningx = binning_esim, - .binningy = {}, - .table = "Threshold", - .x = ".energy", - .y = "", - .v = "ThresholdN", - .err = "ErrThresholdN", - .axisx = "Energy lg(E_{sim}/GeV)", - .axisy = "Counts", - .axisz = "", - .stats = true - }); - - WriteHistogram(connection, { - .name = "SimSpectrumW", - .title = "Weighted Simulated Excess Spectrum", - .dir = "Esim/Simulation/Weighted", - .binningx = binning_esim, - .binningy = {}, - .table = "Threshold", - .x = ".energy", - .y = "", - .v = "Flux", - .err = "ErrFlux", - .axisx = "Energy lg(E_{sim}/GeV)", - .axisy = "dN/dE [cm^{-2} s^{-1} TeV^{-1}]", - .axisz = "", - .stats = true - }); - - 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", - .dir = "Esim", - .binningx = binning_esim, - .binningy = {}, - .table = "Threshold", - .x = ".energy", - .y = "", - .v = "Bias", - .err = "Resolution", - .axisx = "Energy lg(E_{sim}/GeV)", - .axisy = "lg(E_{est}/E_{sim})", - .axisz = "", - .stats = false - }); - - WriteHistogram(connection, { - .name = "Resolution", - .title = "Energy Resolution", - .dir = "Esim", - .binningx = binning_esim, - .binningy = {}, - .table = "Threshold", - .x = ".energy", - .y = "", - .v = "Resolution", - .err = "", - .axisx = "Energy lg(E_{sim}/GeV)", - .axisy = "#sigma(lg(E_{est}/E_{sim}))", - .axisz = "", - .stats = false - }); - - - // ------------------------------------------------------------------- - // ---------------------------- Migration ---------------------------- - // ------------------------------------------------------------------- - - cout << separator("Migration") << '\n'; - - Time start10; - - /* 02:obtain-migration-matrix.sql */ - - // This query gets the analysis results versus Estimated Energy from the combined table - mysqlpp::Query query10(&connection); - query10 << - "CREATE TEMPORARY TABLE Migration ENGINE=Memory AS\n" - "(\n" - " SELECT\n" - " `.energyest`,\n" - " `.energysim`,\n" - " BinningEnergySim.lo AS EsimLo,\n" - " BinningEnergySim.hi AS EsimHi,\n" - " BinningEnergyEst.lo AS EestLo,\n" - " BinningEnergyEst.hi AS EestHi,\n" - " ANY_VALUE(MigrationW) AS MigrationW,\n" - " ANY_VALUE(MigrationN) AS MigrationN\n" - // FIXME: Errors - " FROM\n" - " AnalysisMC\n" - " INNER JOIN\n" - " BinningEnergyEst ON `.energyest`=BinningEnergyEst.bin\n" - " INNER JOIN\n" - " BinningEnergySim ON `.energysim`=BinningEnergySim.bin\n" - " GROUP BY\n" - " `.energyest`, `.energysim`\n" - " ORDER BY\n" - " `.energyest`, `.energysim`\n" - ")"; - - if (print_queries) - PrintQuery(query10.str()); - - qlog << query10 << ";\n" << endl; - if (connection.connected()) - { - cout << query10.execute().info() << endl; - ShowWarnings(connection); - Dump(flog, connection, "Migration"); - - const auto sec10 = Time().UnixTime()-start10.UnixTime(); - cout << "Execution time: " << sec10 << "s\n\n"; - } - - WriteHistogram(connection, { - .name = "MigrationN", - .title = "Energy Migration", - .dir = "", - .binningx = binning_esim, - .binningy = binning_eest, - .table = "Migration", - .x = ".energysim", - .y = ".energyest", - .v = "MigrationN", - .err = "", - .axisx = "Energy lg(E_{sim}/GeV)", - .axisy = "Energy lg(E_{est}/GeV)", - .axisz = "Counts", - .stats = false - }); - - WriteHistogram(connection, { - .name = "MigrationW", - .title = "Energy Migration", - .dir = "", - .binningx = binning_esim, - .binningy = binning_eest, - .table = "Migration", - .x = ".energysim", - .y = ".energyest", - .v = "MigrationW", - .err = "", - .axisx = "Energy lg(E_{sim}/GeV)", - .axisy = "Energy lg(E_{est}/GeV)", - .axisz = "Weighted Counts", - .stats = false - }); - - - - // ------------------------------------------------------------------- - // --------------------------- 11: Summary --------------------------- + } + + // ------------------------------------------------------------------- + // ----------------------------- Summary ----------------------------- + // ------------------------------------------------------------------- cout << separator("Summary") << '\n';