Index: trunk/FACT++/spectrum/analysis.sql =================================================================== --- trunk/FACT++/spectrum/analysis.sql (revision 19882) +++ trunk/FACT++/spectrum/analysis.sql (revision 19882) @@ -0,0 +1,137 @@ +/* ************************************************************************ + + This is the analysis query. It returns two columns for all + signal/background events. + + - `Weight` (Positive for signal (typ. 1), + negative for background (typ. -0.2)) + - `LogEnergyEst` logarithm of estimated energy in GeV + + In additon, all columns provided by the %0:column clause must be + returned. (Note that you must not add a comma behind it) + + %101:files table containing the `FileId`s to analyze. + %102:runinfo table with the run info data + %103:events table with the image parameters + %104:positions table with the source positions in the camera + +*************************************************************************** */ + +WITH Table0 AS +( + SELECT + %100:columns -- this could be removed if we can join events via the same columns (without CorsikaNumResuse) + Weight, + Size, + NumUsedPixels, + NumIslands, + Leakage1, + MeanX, + MeanY, + CosDelta, + SinDelta, + M3Long, + SlopeLong, + Width/Length AS WdivL, + PI()*Width*Length AS Area, + cosa*X - sina*Y AS PX, + cosa*Y + sina*X AS PY + FROM + %101:files + LEFT JOIN + %102:runinfo USING (FileId) + LEFT JOIN + %103:events USING (FileId) -- This could be replaced by a user uploaded temporary table + LEFT JOIN + %104:positions USING (FileId, EvtNumber) + CROSS JOIN + Wobble + WHERE + NumUsedPixels>5.5 + AND + NumIslands<3.5 + AND + Leakage1<0.1 +), + +Table1 AS +( + SELECT + %100:columns + Weight, + Size, CosDelta, SinDelta, M3Long, SlopeLong, Leakage1, WdivL, + MeanX - PX/1.02e0 AS DX, + MeanY - PY/1.02e0 AS DY + FROM + Table0 + WHERE + Area < LOG10(Size)*898e0 - 1535e0 +), + +Table2 AS +( + SELECT + %100:columns + Weight, + Size, CosDelta, SinDelta, DX, DY, M3Long, SlopeLong, Leakage1, WdivL, + SQRT(DX*DX + DY*DY) AS Norm + FROM + Table1 +), + +Table3 AS +( + SELECT + %100:columns + Weight, + Size, M3Long, SlopeLong, Leakage1, WdivL, Norm, + LEAST(GREATEST((CosDelta*DY - SinDelta*DX)/Norm, -1), 1) AS LX, + SIGN(CosDelta*DX + SinDelta*DY) AS Sign + FROM + Table2 +), + +Table5 AS +( + SELECT + %100:columns + Weight, + Size, Leakage1, WdivL, LX, + Norm *0.0117193246260285378e0 AS Dist, + M3Long *Sign*0.0117193246260285378e0 AS M3L, + SlopeLong*Sign/0.0117193246260285378e0 AS Slope + FROM + Table3 +), + +Table6 AS +( + SELECT + %100:columns + Weight, + Size, WdivL, Dist, LX, M3L, Slope, + 1.39252e0 + 0.154247e0*Slope + 1.67972e0*(1-1/(1+4.86232e0*Leakage1)) AS Xi + FROM + Table5 +), + +Table7 AS +( + SELECT + %100:columns + Weight, + Size, Dist, LX, + IF (M3L<-0.07 OR (Dist-0.5e0)*7.2e0-Slope<0, -Xi, Xi) * (1-WdivL) AS Disp + FROM + Table6 +) + +SELECT + %100:columns + Weight, + (Disp*Disp + Dist*Dist - 2*Disp*Dist*SQRT(1-LX*LX)) AS ThetaSq, + %105:estimator AS LogEnergyEst +FROM + Table7 +HAVING + ThetaSq<0.024 Index: trunk/FACT++/spectrum/data.sql =================================================================== --- trunk/FACT++/spectrum/data.sql (revision 19882) +++ trunk/FACT++/spectrum/data.sql (revision 19882) @@ -0,0 +1,25 @@ +WITH Table0 AS +( + SELECT + INTERVAL(LogEnergyEst, %106:bins) AS `.energy`, + COUNT(IF(Weight>0, 1, NULL)) AS `Signal`, + COUNT(IF(Weight<0, 1, NULL))/5 AS `Background` + -- SUM(IF(Weight>0, 1./log10(Size), 0)) AS `Signal`, + -- SUM(IF(Weight<0, 1./log10(Size), 0)) AS `Background` + + -- FIXME: Add excess vs theta + FROM + Excess + GROUP BY + `.energy` + ORDER BY + `.energy` +) + +SELECT + *, + `Signal` - `Background` AS `Excess`, + LiMa(`Signal`, `Background`) AS `Significance`, + ExcErr(`Signal`, `Background`) AS `ErrExcess` +FROM + Table0 Index: trunk/FACT++/spectrum/simulation.sql =================================================================== --- trunk/FACT++/spectrum/simulation.sql (revision 19882) +++ trunk/FACT++/spectrum/simulation.sql (revision 19882) @@ -0,0 +1,87 @@ +WITH Table0 AS +( + SELECT + Weight, + INTERVAL(Zd, %106:theta) AS `.theta`, + INTERVAL(LogEnergyEst, %107:energyest) AS `.energyest`, + INTERVAL(log10(Energy), %108:energysim) AS `.energysim`, + (%109:spectrum)/pow(Energy, SpectralIndex) AS SpectralWeight, -- FIXME: Is this correct for files with different Slopes? + + LogEnergyEst - log10(Energy) AS Residual + FROM + Excess +-- Instead of using %%0:columns, we could join back with the data we need +-- INNER JOIN +-- factmc.EventsMC USING(FileId, EvtNumber, CorsikaNumReuse) +-- INNER JOIN +-- factmc.RunInfoMC USING(FIleId) +), + +Table1 AS +( + SELECT + `.energyest`, + `.energysim`, + + -- Signal, Background, Excess + SUM( IF(Weight>0, ZdWeight*SpectralWeight, 0)) OVER EnergyEst AS `SignalW`, + SUM( IF(Weight<0, ZdWeight*SpectralWeight, 0)) OVER EnergyEst AS `BackgroundW`, + SUM( IF(Weight>0, POW(ErrZdWeight*SpectralWeight,2), 0)) OVER EnergyEst AS `SignalW2`, + SUM( IF(Weight<0, POW(ErrZdWeight*SpectralWeight,2), 0)) OVER EnergyEst AS `BackgroundW2`, + COUNT(IF(Weight>0, 1, NULL)) OVER EnergyEst AS `SignalN`, + COUNT(IF(Weight<0, 1, NULL)) OVER EnergyEst AS `BackgroundN`, + -- SUM( IF(Table9.Weight>0, ThetaHist.Weight*spectrum/Energy, 0)) OVER EnergyEst AS `SignalW`, + -- SUM( IF(Table9.Weight<0, ThetaHist.Weight*spectrum/Energy, 0)) OVER EnergyEst AS `BackgroundW`, + -- SUM( IF(Table9.Weight>0, POW(ThetaHist.Weight*spectrum/Energy,2), 0)) OVER EnergyEst AS `SignalW2`, + -- SUM( IF(Table9.Weight<0, POW(ThetaHist.Weight*spectrum/Energy,2), 0)) OVER EnergyEst AS `BackgroundW2`, + + -- Threshold + SUM( Weight * ZdWeight*SpectralWeight ) OVER EnergySim AS `ThresholdW`, + SUM(POW(Weight * ErrZdWeight*SpectralWeight,2)) OVER EnergySim AS `ThresholdW2`, + SUM( Weight ) OVER EnergySim AS `ThresholdN`, + -- SUM( Table9.Weight/Energy * ThetaHist.Weight*spectrum ) OVER EnergySim AS `ThresholdW`, + -- SUM( POW(Table9.Weight/Energy * ThetaHist.Weight*spectrum,2) ) OVER EnergySim AS `ThresholdW2`, + -- SUM( Table9.Weight/Energy ) OVER EnergySim AS `ThresholdN` + + -- Estimators + SUM(IF(Weight>0, ZdWeight*SpectralWeight, 0)) OVER EnergySim AS SimW, + SUM(IF(Weight>0, Residual* ZdWeight*SpectralWeight, 0)) OVER EnergyEst AS EstSum, + SUM(IF(Weight>0, Residual* ZdWeight*SpectralWeight, 0)) OVER EnergySim AS SimSum, + SUM(IF(Weight>0, POW(Residual,2)*ZdWeight*SpectralWeight, 0)) OVER EnergyEst AS EstSum2, + SUM(IF(Weight>0, POW(Residual,2)*ZdWeight*SpectralWeight, 0)) OVER EnergySim AS SimSum2, + + -- Migration + SUM(Weight * ZdWeight*SpectralWeight) OVER Migration AS `MigrationW`, + SUM(Weight ) OVER Migration AS `MigrationN` + + -- FIXME: Add ExcessN vs theta + -- FIXME: Add ExcessW vs theta + + FROM + Table0 + INNER JOIN + ThetaHist USING(`.theta`) + WINDOW + EnergyEst AS (PARTITION BY `.energyest`), + EnergySim AS (PARTITION BY `.energysim`), + Migration AS (PARTITION BY `.energysim`,`.energyest`) +) + +SELECT DISTINCT + + *, + + `SignalW` - `BackgroundW`/5 AS `ExcessW`, + `SignalN` - `BackgroundN`/5 AS `ExcessN`, + + ExcErr(`SignalW2`, `BackgroundW2`/5) AS `ErrExcessW`, + ExcErr(`SignalN`, `BackgroundN` /5) AS `ErrExcessN`, + + IF(SignalW=0, 0, EstSum / SignalW) AS BiasEst, -- FIMXE: Is NULL better? + IF(SimW =0, 0, SimSum / SimW) AS BiasSim, -- FIMXE: Is NULL better? + + IF(SignalW=0, 0, SQRT(EstSum2/SignalW - POW(EstSum/SignalW, 2))) AS ResolutionEst, + IF(SimW =0, 0, SQRT(SimSum2/SimW - POW(SimSum/SimW, 2))) AS ResolutionSim + +FROM + Table1 Index: trunk/FACT++/src/spectrum.cc =================================================================== --- trunk/FACT++/src/spectrum.cc (revision 19882) +++ trunk/FACT++/src/spectrum.cc (revision 19882) @@ -0,0 +1,1938 @@ +#include +#include +#include + +#include "Database.h" +#include "tools.h" +#include "Time.h" +#include "Configuration.h" + +#ifdef HAVE_HIGHLIGHT +#include "srchilite/sourcehighlight.h" +#include "srchilite/langmap.h" +#endif + +#ifdef HAVE_ROOT +#include "TFile.h" +#include "TH1.h" +#include "TH2.h" +#endif + +using namespace std; +using boost::adaptors::transformed; + +namespace fs = boost::filesystem; + +// ------------------------------- Binning ---------------------------------- + +struct Binning : std::set +{ + bool equidist { false }; + + Binning() { } + Binning(const Binning &m) : std::set(m), equidist(m.equidist) { } + Binning(size_t cnt, double mn, double mx) { set(cnt, mn, mx); } + void set(size_t cnt, double mn, double mx) + { + if (cnt==0) + return; + + if (empty()) + equidist = true; + + const double min = ::min(mn, mx); + const double max = ::max(mn, mx); + + const double stp = (max-min)/cnt; + + for (int i=0; i<=cnt; i++) + emplace(min+i*stp); + } + + void add(double val) + { + emplace(val); + equidist = false; + } + + Binning &operator+=(const vector &v) + { + if (!v.empty()) + { + insert(v.cbegin(), v.cend()); + equidist = false; + } + return *this; + } + + Binning operator+(const vector &v) + { + return Binning(*this) += v; + } + + void add(const double &val) + { + emplace(val); + equidist = false; + } + + string list() const + { + return boost::join(*this | transformed([](double d) { return std::to_string(d); }), ","); + } + + string str() const + { + if (!equidist) + return list(); + + return to_string(size()-1)+":"+to_string(*begin())+","+to_string(*rbegin()); + } + + vector vec() const { return vector(begin(), end()); } + + //double operator[](size_t i) const { return vec().at(i); } +}; + +std::istream &operator>>(std::istream &in, Binning &m) +{ + const istreambuf_iterator eos; + string txt(istreambuf_iterator(in), eos); + + const boost::regex expr( + "([0-9]+)[ /,:;]+" + "([-+]?[0-9]*[.]?[0-9]+([eE][-+]?[0-9]+)?)[ /,:;]+" + "([-+]?[0-9]*[.]?[0-9]+([eE][-+]?[0-9]+)?)" + ); + boost::smatch match; + if (!boost::regex_match(txt, match, expr)) + throw runtime_error("Could not evaluate binning: "+txt); + + m = Binning(stoi(match[1].str()), stof(match[2].str()), stof(match[4].str())); + + return in; +} + +std::ostream &operator<<(std::ostream &out, const Binning &m) +{ + return out << m.str(); +} + +// ---------------------------- Configuration ------------------------------- + +void SetupConfiguration(Configuration &conf) +{ + po::options_description control("Calcsource options"); + control.add_options() + ("uri,u", var() +#if BOOST_VERSION >= 104200 + ->required() +#endif + , "Database link as in\n\tuser:password@server[:port]/database[?compress=0|1].") + ("out,o", var(conf.GetName()), "") + ("dry-run", po_switch(), "Only write the queries to the query.log file. Internally overwrites uri with an empty string.") + ; + + po::options_description binnings("Binnings"); + binnings.add_options() + ("theta", var(Binning(90, 0, 90)), "") + ("theta-bin", vars(), "") + ("esim", var(Binning(15, 2, 5)), "") + ("esim-bin", vars(), "") + ("eest", var(Binning(15, 2, 5)), "") + ("eest-bin", vars(), "") + ; + + po::options_description queries("Queries"); + queries.add_options() + ("analysis", var("analysis.sql"), "") + ("data", var("data.sql"), "") + ("simulation", var("simulation.sql"), "") + ; + + po::options_description preparation("Preparation"); + preparation.add_options() + ("source-key", var(5), "") + ("selector", vars(), "") + ("estimator", var()->required(), "") + ("spectrum", var()->required(), "") + ("env.*", var(), "") + ; + + po::options_description debug("Debug options"); + debug.add_options() + ("print-connection", po_switch(), "Print database connection information") + ("print-queries", po_switch(), "") + ("verbose,v", var(1), "Verbosity (0: quiet, 1: default, 2: more, 3, ...)") + ; + + //po::positional_options_description p; + //p.add("file", 1); // The 1st positional options (n=1) + + conf.AddOptions(control); + conf.AddOptions(binnings); + conf.AddOptions(queries); + conf.AddOptions(preparation); + conf.AddOptions(debug); + //conf.SetArgumentPositions(p); +} + +void PrintUsage() +{ + //78 ............................................................................. + cout << + "spectrum - Calculate a spectrum with classical algorithms\n" + "\n\n" + "Usage: spectrum [-u URI] [options]\n" + "\n" + ; + cout << endl; +} + + + +// ----------------------------- Indentation -------------------------------- + +class sindent : public std::streambuf +{ + std::streambuf *sbuf { nullptr }; + std::ostream *owner { nullptr }; + int lastch { 0 }; // start-of-line + std::string str; + +protected: + virtual int overflow(int ch) + { + if (lastch=='\n' && ch != '\n' ) + sbuf->sputn(str.data(), str.size()); + + return sbuf->sputc(lastch = ch); + } + +public: + explicit sindent(std::streambuf *dest, uint16_t indent=0) + : sbuf(dest), str(indent, ' ') + { + } + + explicit sindent(std::ostream& dest, uint16_t indent=0) + : sbuf(dest.rdbuf()), owner(&dest), str(indent, ' ') + { + owner->rdbuf(this); + } + + virtual ~sindent() + { + if (owner) + owner->rdbuf(sbuf); + } + + void set(uint16_t w) { str.assign(w, ' '); } +}; + +struct indent +{ + uint16_t w; + indent(uint16_t _w=0) : w(_w) { } +}; + +std::ostream &operator<<(std::ostream &out, indent m) +{ + sindent *buf=dynamic_cast(out.rdbuf()); + if (buf) + buf->set(m.w); + return out; +} + +struct separator +{ + string s; + uint16_t n { 0 }; + separator(string _s="") : s(_s) { }; + separator(char c='-', uint16_t _n=78) : s(_n, c), n(_n) { }; +}; + +std::ostream &operator<<(std::ostream &out, separator m) +{ + if (m.n) + out << m.s; + else + { + const uint8_t l = (78-m.s.size())/2-3; + out << string(l<1?1:l, '-') << "={ " << m.s << " }=" << string(l<1?1:l, '-'); + if (m.s.size()%2) + out << '-'; + } + return out; +} + +// ------------------------------- MySQL++ ---------------------------------- + +bool ShowWarnings(Database &connection) +{ + if (!connection.connected()) + return true; + + try + { + const auto resw = + connection.query("SHOW WARNINGS").store(); + + for (size_t i=0; i + struct Value + { + T t { def }; + Value() = default; + Value(const T &_t) : t(_t) { } + operator T() const { return t; } + }; + + string name; + string title; + string dir; + Binning binningx; + Binning binningy; + string table; + string x; + string y; + string v; + string err; + string axisx; + string axisy; + string axisz; + Value stats; +}; + +void WriteHistogram(Database &connection, const Histogram &hist) +{ +#ifdef HAVE_ROOT + if (!connection.connected()) + return; + + mysqlpp::Query query(&connection); + query << "SELECT `%0:x` AS X,%1:y `%2:v` AS V%3:err FROM `%4:table`"; + query.parse(); + + query.template_defaults["table"] = hist.table.c_str(); + + query.template_defaults["x"] = hist.x.c_str(); + query.template_defaults["v"] = hist.v.c_str(); + if (!hist.y.empty()) + query.template_defaults["y"] = (" `"+hist.y+"` AS Y,").c_str(); + if (!hist.err.empty()) + query.template_defaults["err"] = (", `"+hist.err+"` AS E").c_str(); + + // PrintQuery(query.str()); + + const mysqlpp::StoreQueryResult res = query.store(); + + const auto vecx = hist.binningx.vec(); + const auto vecy = hist.binningy.vec(); + + TH1 *h = 0; + + if (hist.y.empty()) + { + h = hist.binningx.equidist ? + new TH1D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.front(), vecx.back()) : + new TH1D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.data()); + } + else + { + if (hist.binningy.equidist) + { + h = hist.binningx.equidist ? + new TH2D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.front(), vecx.back(), vecy.size()-1, vecy.front(), vecy.back()) : + new TH2D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.data(), vecy.size()-1, vecy.front(), vecy.back()); + } + else + { + h = hist.binningx.equidist ? + new TH2D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.front(), vecx.back(), vecy.size()-1, vecy.front(), vecy.back()) : + new TH2D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.data(), vecy.size()-1, vecy.front(), vecy.back()); + } + } + + h->SetXTitle(hist.axisx.c_str()); + h->SetYTitle(hist.axisy.c_str()); + h->SetZTitle(hist.axisz.c_str()); + + h->SetMarkerStyle(kFullDotMedium); + + h->SetStats(hist.stats); + + for (auto ir=res.cbegin(); ir!=res.cend(); ir++) + { + const auto &row = *ir; + + const uint32_t x = row["X"]; + const double v = row["V"]; + if (x>h->GetNbinsX()+1) + continue; + + try + { + const uint32_t y = row["Y"]; + if (y>h->GetNbinsY()+1) + continue; + + h->SetBinContent(x, y, v); + + } + catch (const mysqlpp::BadFieldName &e) + { + h->SetBinContent(x, v); + try + { + h->SetBinError(x, double(row["E"])); + } + catch (const mysqlpp::BadFieldName &e) + { + } + } + } + + gFile->cd(); + TDirectory *dir = gFile->GetDirectory(hist.dir.c_str()); + if (dir) + dir->cd(); + else + { + gFile->mkdir(hist.dir.c_str()); + gFile->cd(hist.dir.c_str()); + } + h->Write(); + delete h; +#endif +} + +// -------------------------------- Resources ------------------------------- + +#define RESOURCE(TYPE,RC) ([]() { \ + extern const char _binary_##RC##_start[]; \ + extern const char _binary_##RC##_end[]; \ + return TYPE(_binary_##RC##_start, _binary_##RC##_end); \ +})() + +string CreateResource(Configuration &conf, const string id, const string def, const string resource) +{ + string file = conf.Get(id); + + if (file!=def) + { + file = conf.GetPrefixedString(id); + cout << "Using " << file << "." << endl; + return file; + } + + file = conf.GetPrefixedString(id); + + cout << "Searching " << file << "... "; + + if (fs::exists(file)) + { + cout << "found." << endl; + return file; + } + + cout << "creating!" << endl; + + ofstream(file) << resource; + + return file; +} + +// ================================== MAIN ================================== + + +int main(int argc, const char* argv[]) +{ + + Time start; + + Configuration conf(argv[0]); + conf.SetPrintUsage(PrintUsage); + SetupConfiguration(conf); + + if (!conf.DoParse(argc, argv)) + return 127; + + // ----------------------------- Evaluate options -------------------------- + + const string uri = conf.Get("dry-run") ? "" : conf.Get("uri"); + const string out = conf.Get("out"); + const uint16_t verbose = conf.Get("verbose"); + + 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"); + + cout << "B\"theta\": " << binning_theta.str() << endl; + cout << "B\"esim\": " << binning_esim.str() << endl; + cout << "B\"eest\": " << binning_eest.str() << endl; + + const uint16_t source_key = conf.Get("source-key"); + const string where = boost::join(conf.Vec("selector"), " AND\n "); + const string estimator = conf.Get("estimator"); + const string spectrum = conf.Get("spectrum"); + const auto env = conf.GetOptions("env."); + + cout << "\n"; + const string analysis_sql = CreateResource(conf, "analysis", "analysis.sql", RESOURCE(std::string, spectrum_analysis_sql)); + const string data_sql = CreateResource(conf, "data", "data.sql", RESOURCE(std::string, spectrum_data_sql)); + const string simulation_sql = CreateResource(conf, "simulation", "simulation.sql", RESOURCE(std::string, spectrum_simulation_sql)); + cout << endl; + + const string strzd = binning_theta.list(); + const string stre = binning_eest.list(); + const string strth = binning_esim.list(); + + // ------------------------------------------------------------------------- + // Checking for database connection + + Database connection(uri); // Keep alive while fetching rows + + if (!uri.empty()) + { + if (verbose>0) + { + cout << "Connecting to database...\n"; + cout << "Client Version: " << mysqlpp::Connection().client_version() << endl; + } + + try + { + connection.connected(); + } + catch (const exception &e) + { + cerr << "SQL connection failed: " << e.what() << endl; + return 1; + } + + if (verbose>0) + { + cout << "Server Version: " << connection.server_version() << endl; + cout << "Connected to " << connection.ipc_info() << endl; + } + + if (print_connection) + { + try + { + const auto res1 = connection.query("SHOW STATUS LIKE 'Compression'").store(); + cout << "Compression of database connection is " << string(res1[0][1]) << endl; + + const auto res2 = connection.query("SHOW STATUS LIKE 'Ssl_cipher'").store(); + cout << "Connection to databases is " << (string(res2[0][1]).empty()?"UNENCRYPTED":"ENCRYPTED ("+string(res2[0][1])+")") << endl; + } + catch (const exception &e) + { + cerr << "\nSHOW STATUS LIKE 'Compression'\n\n"; + cerr << "SQL query failed:\n" << e.what() << endl; + return 9; + } + } + } + + // ------------------------------------------------------------------------- + // Create log streams + + ofstream qlog(out+".query.sql"); + ofstream flog(connection.connected() ? out+".dump.sql" : ""); + ofstream mlog(connection.connected() ? out+".C" : ""); + + cout << "\n"; + cout << "Queries will be logged to " << out << ".query.sql\n"; + if (connection.connected()) + { + cout << "Tables will be dumped to " << out << ".dump.sql\n"; + cout << "ROOT macro will be written to " << out << ".C\n"; + } + +#ifdef HAVE_ROOT + TFile root(connection.connected() ? (out+".hist.root").c_str() : "", "RECREATE"); + if (connection.connected()) + { + if (root.IsZombie()) + return 10; + cout << "Histograms will be written to " << out << ".hist.root\n"; + } +#endif + + cout << endl; + + // FIMXE: Implement SYNTAX check on spectrum, estimator and selector + + // ------------------------------------------------------------------- + // --------------------------- 0: Binnings --------------------------- + + cout << separator("Binnings") << '\n'; + + CreateBinning(connection, qlog, binning_theta, "Theta"); + CreateBinning(connection, qlog, binning_eest, "EnergyEst"); + CreateBinning(connection, qlog, binning_esim, "EnergySim"); + + Dump(flog, connection, "BinningTheta"); + Dump(flog, connection, "BinningEnergyEst"); + Dump(flog, connection, "BinningEnergySim"); + + // ------------------------------------------------------------------- + // --------------------------- 1: DataFiles -------------------------- + cout << separator("DataFiles") << '\n'; + + Time start1; + + /* 01:get-data-files.sql */ + mysqlpp::Query query1(&connection); + query1 << + "CREATE TEMPORARY TABLE DataFiles\n" + "(\n" + " FileId INT UNSIGNED NOT NULL,\n" + " PRIMARY KEY (FileId) USING HASH\n" + ") ENGINE=Memory\n" + "AS\n" + "(\n" + " SELECT\n" + " FileId\n" + " FROM\n" + " factdata.RunInfo\n" + " WHERE\n" + " fRunTypeKEY=1 AND fSourceKEY=%100:source AND\n" + " %101:where\n" + " ORDER BY\n" + " FileId\n" // In order: faster + ")"; + + query1.parse(); + for (auto it=env.cbegin(); it!=env.cend(); it++) + query1.template_defaults[it->first.c_str()] = it->second.c_str(); + + query1.template_defaults["source"] = to_string(source_key).c_str(); + query1.template_defaults["where"] = where.c_str(); + + if (print_queries) + PrintQuery(query1.str()); + + qlog << query1 << ";\n" << endl; + if (connection.connected()) + { + cout << query1.execute().info() << endl; + ShowWarnings(connection); + Dump(flog, connection, "DataFiles"); + + const auto sec1 = Time().UnixTime()-start1.UnixTime(); + cout << "Execution time: " << sec1 << "s\n" << endl; + } + + // FIXME: Setup Zd binning depending on Data + + // ------------------------------------------------------------------- + // ------------------------ 2: ObservationTime ----------------------- + + cout << separator("ObservationTime") << '\n'; + + Time start2; + + /* 02:get-observation-time.sql */ + mysqlpp::Query query2(&connection); + query2 << + "CREATE TEMPORARY TABLE ObservationTime\n" + "(\n" + " `.theta` SMALLINT UNSIGNED NOT NULL,\n" + " OnTime FLOAT NOT NULL,\n" + " PRIMARY KEY (`.theta`) USING HASH\n" + ") ENGINE=Memory\n" + "AS\n" + "(\n" + " SELECT\n" + " INTERVAL(fZenithDistanceMean, %100:bins) AS `.theta`,\n" + " SUM(TIME_TO_SEC(TIMEDIFF(fRunStop,fRunStart))*fEffectiveOn) AS OnTime\n" + " FROM\n" + " DataFiles\n" + " LEFT JOIN \n" + " factdata.RunInfo USING (FileId)\n" + " GROUP BY\n" + " `.theta`\n" + " ORDER BY\n" + " `.theta`\n" + ")"; + + query2.parse(); + for (auto it=env.cbegin(); it!=env.cend(); it++) + query2.template_defaults[it->first.c_str()] = it->second.c_str(); + + query2.template_defaults["bins"] = strzd.c_str(); + + if (print_queries) + PrintQuery(query2.str()); + + qlog << query2 << ";\n" << endl; + if (connection.connected()) + { + cout << query2.execute().info() << endl; + ShowWarnings(connection); + Dump(flog, connection, "ObservationTime"); + + const auto sec2 = Time().UnixTime()-start2.UnixTime(); + cout << "Execution time: " << sec2 << "s\n\n"; + } + + /* + + SELECT + min(first_number) as first_number, + max(last_number) as last_number + from + ( + select + first_number, + last_number + sum(grp) over(order by first_number) as grp + from + ( + select + bin AS first_number, + LEAD(bin) AS last_number + IF(first_number <> lag(bin, 2) over (order by bin) + 1, 1, 0) as grp + from + t1 + ) + ) + group by grp + order by 1 + */ + + /* +WITH Hallo AS (WITH MyTest AS (WITH MyTable AS (SELECT bin FROM Test +WHERE cnt>0 GROUP BY bin ORDER BY bin) SELECT bin AS first_number, +LAG(bin) OVER(ORDER BY bin) AS prev, LEAD(bin) OVER(ORDER BY bin) AS +last_number, IF(bin <> (LAG(bin) OVER(ORDER BY bin)+1), 1, 0) AS +grp_prev, IF(bin <> (LEAD(bin) OVER(ORDER BY bin)-1), 1, 0) as grp_nxt +FROM MyTable ) SELECT first_number, last_number, sum(grp_nxt) +over(order by first_number) as grp FROM MyTest) SELECT * FROM Hallo + + + WITH MyTable AS +(SELECT bin FROM Test WHERE cnt>0 GROUP BY bin ORDER BY bin) +SELECT bin, LAG(bin) OVER(ORDER BY bin) AS prev, LEAD(bin) OVER(ORDER BY bin) AS nxt, +IF(bin = (LAG(bin) OVER(ORDER BY bin)+1), 1, 0) AS grp_prev, IF(bin = (LEAD(bin) OVER(ORDER BY bin)-1), 1, 0) as grp_nxt FROM MyTable + +"-1","0" +"0","0" +"2","3" 2-4 +"3","5" +"4","6" +"5","0" +"6","0" +"7","7" 7-7 +"8","0" +"9","9" 9-12 +"10","10" +"11","11" +"12","12" +"13","0" +"14","0" +"15","0" +"16","16" 16-19 +"17","17" +"18","18" +"19","19" + */ + + // ------------------------------------------------------------------- + // ------------------------ 3: MonteCarloFiles ----------------------- + + cout << separator("MonteCarloFiles") << '\n'; + + Time start3; + + /* 02:get-monte-carlo-files.sql */ + mysqlpp::Query query3(&connection); + query3 << + "CREATE TEMPORARY TABLE MonteCarloFiles\n" + "(\n" + " FileId INT UNSIGNED NOT NULL,\n" + " PRIMARY KEY (FileId) USING HASH\n" + ") ENGINE=Memory\n" + "AS\n" + "(\n" + " SELECT\n" + " FileId\n" + " FROM\n" + " ObservationTime\n" + " LEFT JOIN\n" + " BinningTheta ON `.theta`=bin\n" + " LEFT JOIN\n" + " factmc.RunInfoMC\n" + " ON\n" + //" ThetaMin BETWEEN lo AND hi OR ThetaMax BETWEEN lo AND hi\n" // Includes BOTH edges + " (ThetaMin>=lo AND ThetaMinlo AND ThetaMax<=hi)\n" + " WHERE\n" + " PartId=1 AND\n" + " FileId%%2=0\n" + " ORDER BY\n" + " FileId\n" // In order: faster + ")"; +/* + " SELECT\n" + " FileId\n" + " FROM\n" + " factmc.RunInfoMC\n" + " WHERE\n" + " PartId=1 AND\n" + " FileId%%2=0 AND\n" + " (%0:range)\n" + ")"; +*/ + + query3.parse(); + for (auto it=env.cbegin(); it!=env.cend(); it++) + query3.template_defaults[it->first.c_str()] = it->second.c_str(); + + if (print_queries) + PrintQuery(query3.str()); + + qlog << query3 << ";\n" << endl; + if (connection.connected()) + { + cout << query3.execute().info() << endl; + ShowWarnings(connection); + Dump(flog, connection, "MonteCarloFiles"); + + const auto sec3 = Time().UnixTime()-start3.UnixTime(); + cout << "Execution time: " << sec3 << "s\n\n"; + } + + // ------------------------------------------------------------------- + // ----------------------- 3b: Monte Carlo Area ---------------------- + + cout << separator("MonteCarloArea") << '\n'; + + Time start3b; + + /* 02:get-monte-carlo-files.sql */ + mysqlpp::Query query3b(&connection); + query3b << + "CREATE TEMPORARY TABLE MonteCarloArea ENGINE=Memory\n" + "AS\n" + "(\n" + " SELECT\n" + " MIN(`CSCAT[1]`) AS MinImpactLo,\n" + " MAX(`CSCAT[1]`) AS MaxImpactLo,\n" + " MIN(`CSCAT[2]`) AS MinImpactHi,\n" + " MAX(`CSCAT[2]`) AS MaxImpactHi\n" + " FROM\n" + " MonteCarloFiles\n" + " LEFT JOIN\n" + " factmc.CorsikaSetup ON FileId=RUNNR\n" + " GROUP BY\n" + " `CSCAT[1]`, `CSCAT[2]`\n" + " ORDER BY\n" + " MaxImpactHi\n" + ")"; + + query3b.parse(); + for (auto it=env.cbegin(); it!=env.cend(); it++) + query3b.template_defaults[it->first.c_str()] = it->second.c_str(); + + if (print_queries) + PrintQuery(query3b.str()); + + qlog << query3b << ";\n" << endl; + if (connection.connected()) + { + cout << query3b.execute().info() << endl; + ShowWarnings(connection); + if (Dump(flog, connection, "MonteCarloArea")!=1) + { + cerr << "Impact range inconsistent!" << endl; + return 1; + } + + const auto sec3b = Time().UnixTime()-start3b.UnixTime(); + cout << "Execution time: " << sec3b << "s\n\n"; + } + + + // ------------------------------------------------------------------- + // -------------------------- 4: ThetaHist --------------------------- + + cout << separator("ThetaHist") << '\n'; + + Time start4; + + /* 02:get-theta-distribution.sql */ + mysqlpp::Query query4(&connection); + query4 << + "CREATE TEMPORARY TABLE ThetaHist\n" + "(\n" + " `.theta` SMALLINT UNSIGNED NOT NULL,\n" + " lo DOUBLE NOT NULL COMMENT 'Lower edge of zenith distance bin in degree',\n" + " hi DOUBLE NOT NULL COMMENT 'Upper edge of zenith distance bin in degree',\n" + " CountN INT UNSIGNED NOT NULL,\n" + " ErrCountN DOUBLE NOT NULL,\n" + " OnTime FLOAT NOT NULL,\n" + " ZdWeight DOUBLE NOT NULL COMMENT 'tau(delta theta)',\n" + " ErrZdWeight DOUBLE NOT NULL COMMENT 'sigma(tau)',\n" + " PRIMARY KEY (`.theta`) USING HASH\n" + ") ENGINE=Memory\n" + "AS\n" + "(\n" + " WITH EventCount AS\n" + " (\n" + " SELECT\n" + " INTERVAL(DEGREES(Theta), %100:bins) AS `.theta`,\n" + " COUNT(*) AS CountN\n" + " FROM\n" + " MonteCarloFiles\n" + " LEFT JOIN\n" + " factmc.OriginalMC USING(FileId)\n" + " GROUP BY\n" + " `.theta`\n" + " )\n" + " SELECT\n" + " `.theta`, lo, hi,\n" + " CountN,\n" + " SQRT(CountN) AS ErrCountN,\n" + " OnTime,\n" + " OnTime/CountN AS ZdWeight,\n" + " (OnTime/CountN)*SQRT(POW(1/300, 2) + 1/CountN) AS ErrZdWeight\n" + " FROM\n" + " ObservationTime\n" + " LEFT JOIN\n" + " EventCount USING(`.theta`)\n" + " LEFT JOIN\n" + " BinningTheta ON `.theta`=bin\n" + " ORDER BY\n" + " `.theta`\n" + ")"; + + query4.parse(); + for (auto it=env.cbegin(); it!=env.cend(); it++) + query4.template_defaults[it->first.c_str()] = it->second.c_str(); + + query4.template_defaults["bins"] = strzd.c_str(); + + if (print_queries) + PrintQuery(query4.str()); + + qlog << query4 << ";\n" << endl; + if (connection.connected()) + { + cout << query4.execute().info() << endl; + ShowWarnings(connection); + Dump(flog, connection, "ThetaHist"); + + const auto sec4 = Time().UnixTime()-start4.UnixTime(); + cout << "Execution time: " << sec4 << "s\n"; + + + const mysqlpp::StoreQueryResult res4 = connection.query("SELECT * FROM ThetaHist").store(); + + cout << " Bin MC Counts OnTime ZdWeight\n"; + const auto bins = binning_theta.vec(); + for (auto ir=res4.cbegin(); ir!=res4.cend(); ir++) + { + const mysqlpp::Row &row = *ir; + + const uint32_t bin = row[".theta"]; + cout << setw(5) << bins[bin] << ": " << setw(10) << row["CountN"] << " " << setw(10) << row["OnTime"] << " " << setw(10) << row["ZdWeight"] << '\n'; + } + cout << endl; + } + + WriteHistogram(connection, { + .dir = "Zd", + .name = "OnTime", + .title = "Effective on time", + .binningx = binning_theta, + .table = "ThetaHist", + .x = ".theta", + .v = "OnTime", + .axisx = "Zenith Distance #theta [#circ]", + .axisy = "Eff. on time [s]" + }); + + WriteHistogram(connection, { + .dir = "Zd", + .name = "CountN", + .title = "Simulated Zenith Distance", + .binningx = binning_theta, + .table = "ThetaHist", + .x = ".theta", + .v = "CountN", + .err = "ErrCountN", + .axisx = "Zenith Distance #theta [#circ]", + .axisy = "Counts" + }); + + WriteHistogram(connection, { + .dir = "Zd", + .name = "ZdWeight", + .title = "Zenith Distance Weight", + .binningx = binning_theta, + .table = "ThetaHist", + .x = ".theta", + .v = "ZdWeight", + .err = "ErrZdWeight", + .axisx = "Zenith Distance #theta [#circ]", + .axisy = "Weight [s]" + }); + + + // ------------------------------------------------------------------- + // ------------------------- 5: AnalysisData ------------------------- + + cout << separator("AnalysisData") << '\n'; + + Time start5; + + /* 02:analyze-data.sql */ + mysqlpp::Query query5(&connection); + sindent indent5(query5); + query5 << + "CREATE TEMPORARY TABLE AnalysisData\n" + "(\n" + " `.energy` SMALLINT UNSIGNED NOT NULL,\n" + " `Signal` DOUBLE NOT NULL,\n" + " `Background` DOUBLE NOT NULL,\n" + " `Excess` DOUBLE NOT NULL,\n" + " `Significance` DOUBLE NOT NULL,\n" + " `ErrExcess` DOUBLE NOT NULL,\n" + " PRIMARY KEY (`.energy`) USING HASH\n" + ") ENGINE=Memory\n" + "AS\n" + "(\n" + " WITH Excess AS\n" + " (\n" << indent(6) + << ifstream(analysis_sql).rdbuf() << indent(0) << + " ),\n" << + " Result AS\n" + " (\n" << indent(6) + << ifstream(data_sql).rdbuf() << indent(0) << // Must end with EOL and not in the middle of a comment + " )\n" + " SELECT * FROM Result\n" + ")"; + + query5.parse(); + for (auto it=env.cbegin(); it!=env.cend(); it++) + query5.template_defaults[it->first.c_str()] = it->second.c_str(); + + //query5.template_defaults["columns"] = "FileId, EvtNumber,"; + query5.template_defaults["columns"] = ""; + query5.template_defaults["files"] = "DataFiles"; + query5.template_defaults["runinfo"] = "factdata.RunInfo"; + query5.template_defaults["events"] = "factdata.Images"; + query5.template_defaults["positions"] = "factdata.Position"; + + query5.template_defaults["bins"] = stre.c_str(); + query5.template_defaults["estimator"] = estimator.c_str(); + + if (print_queries) + PrintQuery(query5.str()); + + qlog << query5 << ";\n" << endl; + if (connection.connected()) + { + cout << query5.execute().info() << endl; + ShowWarnings(connection); + Dump(flog, connection, "AnalysisData"); + + const auto sec5 = Time().UnixTime()-start5.UnixTime(); + cout << "Execution time: " << sec5 << "s\n"; + + + const mysqlpp::StoreQueryResult res5 = connection.query("SELECT * FROM AnalysisData").store(); + + cout << " Bin Signal Background Excess Significance Error" << endl; + for (auto row=res5.cbegin(); row!=res5.cend(); row++) + { + for (auto it=row->begin(); it!=row->end(); it++) + cout << setw(10) << *it << " "; + cout << '\n'; + } + cout << endl; + } + + + // ------------------------------------------------------------------- + // --------------------------- 6: ResultMC --------------------------- + + cout << separator("ResultMC") << '\n'; + + Time start6; + + /* 02:analyze-simulation.sql */ + + // This table combines the analysis results vs. Binning in Estimated Energy and Simulated Energy + mysqlpp::Query query6(&connection); + sindent indent6(query6); + query6 << + "CREATE TEMPORARY TABLE AnalysisMC\n" + "(\n" + " `.energyest` SMALLINT UNSIGNED NOT NULL,\n" + " `.energysim` SMALLINT UNSIGNED NOT NULL,\n" + " SignalW DOUBLE NOT NULL,\n" + " BackgroundW DOUBLE NOT NULL,\n" + " ThresholdW DOUBLE NOT NULL,\n" + " SignalN DOUBLE NOT NULL,\n" + " BackgroundN DOUBLE NOT NULL,\n" + " ThresholdN DOUBLE NOT NULL,\n" + " ExcessW DOUBLE NOT NULL,\n" + " ExcessN DOUBLE NOT NULL,\n" + " ErrExcessW DOUBLE NOT NULL,\n" + " ErrExcessN DOUBLE NOT NULL,\n" + " BiasEst DOUBLE NOT NULL,\n" + " BiasSim DOUBLE NOT NULL,\n" + " ResolutionEst DOUBLE,\n" + " ResolutionSim DOUBLE,\n" + " INDEX (`.energyest`) USING HASH,\n" + " INDEX (`.energysim`) USING HASH\n" + ") ENGINE=Memory\n" + "AS\n" + "(\n" + " WITH Excess AS\n" + " (\n" << indent(6) + << ifstream(analysis_sql).rdbuf() << indent(0) << + " ),\n" << + " Result AS\n" + " (\n" << indent(6) + << ifstream(simulation_sql).rdbuf() << indent(0) << // Must end with EOL and not in the middle of a comment + " )\n" + " SELECT * FROM Result\n" + ")"; + + query6.parse(); + for (auto it=env.cbegin(); it!=env.cend(); it++) + query6.template_defaults[it->first.c_str()] = it->second.c_str(); + + //query6.template_defaults["columns"] = "FileId, EvtNumber, CorsikaNumReuse,"; + query6.template_defaults["columns"] = "Zd, Energy, SpectralIndex,"; + query6.template_defaults["files"] = "MonteCarloFiles"; + query6.template_defaults["runinfo"] = "factmc.RunInfoMC"; + query6.template_defaults["events"] = "factmc.EventsMC"; + query6.template_defaults["positions"] = "factmc.PositionMC"; + + query6.template_defaults["energyest"] = stre.c_str(); + query6.template_defaults["energysim"] = strth.c_str(); + query6.template_defaults["theta"] = strzd.c_str(); + query6.template_defaults["spectrum"] = spectrum.c_str(); + query6.template_defaults["estimator"] = estimator.c_str(); + + if (print_queries) + PrintQuery(query6.str()); + + qlog << query6 << ";\n" << endl; + if (connection.connected()) + { + cout << query6.execute().info() << endl; + ShowWarnings(connection); + Dump(flog, connection, "AnalysisMC"); + + const auto sec6 = Time().UnixTime()-start6.UnixTime(); + cout << "Execution time: " << sec6 << "s\n\n"; + } + + + // ------------------------------------------------------------------- + // ----------------------- 7: SimulatedSpectrum ---------------------- + + cout << separator("SimulatedSpectrum") << '\n'; + + // [lo^(1+gammaa) - hi^(1+gamma)] / (1+gamma) + + Time start7; + /* 02:get-corsika-events.sql */ + + // FIXME: Theta weights? + // FIXME: energysim binning + mysqlpp::Query query7(&connection); + query7 << + "CREATE TEMPORARY TABLE SimulatedSpectrum\n" + "(\n" + " `.energy` SMALLINT UNSIGNED NOT NULL COMMENT 'Bin Index [MC Energy]',\n" + " CountN DOUBLE NOT NULL,\n" // COMMENT 'Number of events',\n" + " CountW DOUBLE NOT NULL,\n" // COMMENT 'Sum of weights, reweighted sum',\n" + " CountW2 DOUBLE NOT NULL,\n" // COMMENT 'Sum of square of weights'\n" + " PRIMARY KEY (`.energy`) USING HASH\n" + ") ENGINE=Memory\n" + "AS\n" + "(\n" + " SELECT\n" + " INTERVAL(LOG10(Energy), %100:energyest) AS `.energy`,\n" + " COUNT(*) AS CountN,\n" + " SUM( (%101:spectrum)/pow(Energy, SpectralIndex) ) AS CountW,\n" // Contents is: CountW + " SUM(POW((%101:spectrum)/pow(Energy, SpectralIndex),2)) AS CountW2\n" // Error is: SQRT(CountW2) + " FROM\n" + " MonteCarloFiles\n" + " LEFT JOIN\n" + " factmc.RunInfoMC USING (FileId)\n" + " LEFT JOIN\n" + " factmc.OriginalMC USING (FileId)\n" + " GROUP BY\n" + " `.energy`\n" + " ORDER BY\n" + " `.energy`\n" + ")"; + + query7.parse(); + for (auto it=env.cbegin(); it!=env.cend(); it++) + query7.template_defaults[it->first.c_str()] = it->second.c_str(); + + //query7.template_defaults["list"] = listmc.c_str(); + query7.template_defaults["energyest"] = stre.c_str(); + query7.template_defaults["spectrum"] = spectrum.c_str(); + + if (print_queries) + PrintQuery(query7.str()); + + qlog << query7 << ";\n" << endl; + if (connection.connected()) + { + cout << query7.execute().info() << endl; + ShowWarnings(connection); + Dump(flog, connection, "SimulatedSpectrum"); + + const auto sec7 = Time().UnixTime()-start7.UnixTime(); + cout << "Execution time: " << sec7 << "s\n"; + + + const mysqlpp::StoreQueryResult res7 = connection.query("SELECT * FROM SimulatedSpectrum").store(); + + //cout << "Received rows: " << res7.num_rows() << '\n' << endl; + + // 123456789|123456789|123456789|123456789|123456789|123456789|123456789 + cout << " Bin CountW CountW2" << endl; + const auto bins = binning_esim.vec(); + for (auto ir=res7.cbegin(); ir!=res7.cend(); ir++) + { + const mysqlpp::Row &row = *ir; + + const uint32_t bin = row[".energy"]; + cout << setw(5) << bins[bin] << ": " << setw(10) << row["CountW"] << " " << setw(10) << row["CountW2"] << '\n'; + } + cout << endl; + } + + + // ------------------------------------------------------------------- + // --------------------------- 8: Spectrum --------------------------- + + cout << separator("Spectrum") << '\n'; + + Time start8; + + /* 02:calculate-spectrum.sql */ + + mysqlpp::Query query8(&connection); + query8 << + // SELECT SUM(CountN) AS CountN, SUM(OnTime) AS OnTime FROM ThetaHist + "CREATE TEMPORARY TABLE Spectrum\n" + "(\n" + " `.energy` SMALLINT UNSIGNED NOT NULL COMMENT 'Bin Index [Energy]' PRIMARY KEY,\n" + " lo DOUBLE NOT NULL COMMENT 'Lower edge of energy bin in lg(E/GeV)',\n" + " hi DOUBLE NOT NULL COMMENT 'Upper edge of energy bin in lg(E/GeV)',\n" + " `Signal` DOUBLE NOT NULL COMMENT 'Number of signal events',\n" + " `Background` DOUBLE NOT NULL COMMENT 'Average number of background events',\n" + " `Excess` DOUBLE NOT NULL COMMENT 'Number of excess events',\n" + " ErrSignal DOUBLE NOT NULL COMMENT 'Poisson error on number of signal events',\n" + " ErrBackground DOUBLE NOT NULL COMMENT 'Poisson error on number of background events',\n" + " `ErrExcess` DOUBLE NOT NULL COMMENT 'Error of excess events',\n" + " `Significance` DOUBLE NOT NULL COMMENT 'Li/Ma sigficance',\n" + " `ExcessN` DOUBLE NOT NULL COMMENT 'Number of excess events in simulated data',\n" + " `ExcessW` DOUBLE NOT NULL COMMENT 'Weighted number of excess events in simulated data',\n" + " `ErrExcessN` DOUBLE NOT NULL COMMENT 'Error or number of excess events in simulated data',\n" + " `ErrExcessW` DOUBLE NOT NULL COMMENT 'Error of weighted number of excess events in simulated data',\n" + " SignalW DOUBLE NOT NULL COMMENT 'Weighted number of signal events in simulated data',\n" + " BackgroundW DOUBLE NOT NULL COMMENT 'Weighted number of background events in simulated data',\n" + " ErrSignalW DOUBLE NOT NULL COMMENT 'Error of weighted number of signal events in simulated data',\n" + " ErrBackgroundW DOUBLE NOT NULL COMMENT 'Error of weighted number of background events in simulated data',\n" + " Flux DOUBLE NOT NULL COMMENT 'dN/dA/dt [cm^-2 s-^1 TeV^-1]',\n" + " ErrFlux DOUBLE NOT NULL COMMENT 'dN/dA/dt [cm^-2 s-^1 TeV^-1]',\n" + " Bias DOUBLE NOT NULL COMMENT 'Energy Bias, average residual in lg(E)',\n" + " Resolution DOUBLE NOT NULL COMMENT 'Energy resolution, standard divation of residual in lg(E)',\n" + " EfficiencyN DOUBLE NOT NULL COMMENT 'Simulated cut efficiency (weighted)',\n" + " EfficiencyW DOUBLE NOT NULL COMMENT 'Simulated cut efficiency (unweighted)',\n" + " ErrEfficiencyN DOUBLE NOT NULL COMMENT 'Error of simulated cut efficiency (weighted)',\n" + " ErrEfficiencyW DOUBLE NOT NULL COMMENT 'Error of simulated cut efficiency (unweighted)'\n" + ") ENGINE=Memory\n" + "AS\n" + "(\n" + " WITH ThetaSums AS\n" + " (\n" + " SELECT\n" + " SUM(CountN) AS CountSim,\n" + " SUM(OnTime) AS ObsTime\n" + " FROM\n" + " ThetaHist\n" + " ),\n" + " Area AS\n" + " (\n" + " SELECT\n" + " POW(MinImpactHi,2)*PI() AS Area\n" + " FROM\n" + " MonteCarloArea\n" + " ),\n" + " ResultMC AS\n" // Group AnalysisMC by EnergyEst Binning + " (\n" + " SELECT\n" + " `.energyest` AS `.energy`,\n" + " ANY_VALUE(SignalW) AS SignalW,\n" + " ANY_VALUE(SignalW2) AS SignalW2,\n" + " ANY_VALUE(BackgroundW) AS BackgroundW,\n" + " ANY_VALUE(BackgroundW2) AS BackgroundW2,\n" + " ANY_VALUE(SignalN) AS SignalN,\n" + " ANY_VALUE(BackgroundN) AS BackgroundN,\n" + " ANY_VALUE(ExcessW) AS ExcessW,\n" + " ANY_VALUE(ExcessN) AS ExcessN,\n" + " ANY_VALUE(ErrExcessW) AS ErrExcessW,\n" + " ANY_VALUE(ErrExcessN) AS ErrExcessN,\n" + " ANY_VALUE(BiasEst) AS Bias,\n" + " ANY_VALUE(ResolutionEst) AS Resolution\n" + " FROM\n" + " AnalysisMC\n" + " GROUP BY\n" + " `.energy`\n" + " ORDER BY\n" + " `.energy`\n" + " )\n" + " SELECT\n" + " `.energy`, lo, hi,\n" // Scale for Theta-Weights + " `Signal`, `Background`/5 AS `Background`, `Excess`, `ErrExcess`, `Significance`,\n" + " SQRT(`Signal`) AS ErrSignal,\n" + " SQRT(`SignalW2`) AS ErrSignalW,\n" + " SQRT(`Background`)/5 AS ErrBackground,\n" + " SQRT(`BackgroundW2`)/5 AS ErrBackgroundW,\n" + " ExcessN, ExcessW, ErrExcessN, ErrExcessW, SignalW, BackgroundW,\n" + " AnalysisData.Excess/ResultMC.ExcessW*SimulatedSpectrum.CountW * 1000/(POW(10,hi)-POW(10,lo)) /Area/ObsTime / CountSim*ObsTime AS Flux,\n" + " AnalysisData.Excess/ResultMC.ExcessW*SimulatedSpectrum.CountW * 1000/(POW(10,hi)-POW(10,lo)) /Area/ObsTime / CountSim*ObsTime\n" + " * SQRT(\n" + " + POW(AnalysisData.ErrExcess / AnalysisData.Excess, 2)\n" + " + POW(ResultMC.ErrExcessW / ResultMC.ExcessW, 2)\n" + " + SimulatedSpectrum.CountW2 / POW(SimulatedSpectrum.CountW,2)\n" + " ) AS ErrFlux,\n" + " Bias,\n" + " Resolution,\n" + " ResultMC.ExcessW/SimulatedSpectrum.CountW * CountSim/ObsTime AS EfficiencyW,\n" + " ResultMC.ExcessN/SimulatedSpectrum.CountN AS EfficiencyN,\n" + " ( POW(ResultMC.ErrExcessW/ResultMC.ExcessW, 2) + POW(SQRT(SimulatedSpectrum.CountW2)/SimulatedSpectrum.CountW, 2) )\n" + " * POW(ResultMC.ExcessW/SimulatedSpectrum.CountW, 2) * CountSim/ObsTime AS ErrEfficiencyW,\n" + " ( POW(ResultMC.ErrExcessN, 2) + POW(ResultMC.ExcessN, 2)/SimulatedSpectrum.CountN)/POW(SimulatedSpectrum.CountN, 2) AS ErrEfficiencyN\n" + " FROM\n" + " AnalysisData\n" + " INNER JOIN\n" + " ResultMC USING(`.energy`)\n" + " INNER JOIN\n" + " SimulatedSpectrum USING(`.energy`)\n" + " INNER JOIN\n" + " BinningEnergyEst ON `.energy`=bin\n" + " CROSS JOIN\n" + " ThetaSums, Area\n" + " WHERE\n" + " AnalysisData.Excess>0\n" + " ORDER BY\n" + " `.energy`\n" + ")" + ; + + // [ Sa^2/a^2 + Sb^2/b^2 ] * a^2/b^2 + // [ (sc^2)/c^2+(sb^2)/b^2+(sa^2)/a^2 ] * a^2*b^2/c^2 + + query8.parse(); + for (auto it=env.cbegin(); it!=env.cend(); it++) + query8.template_defaults[it->first.c_str()] = it->second.c_str(); + + //query8.template_defaults["area"] = area; + //query8.template_defaults["ontime"] = resX[0]["OnTime"].data(); + //query8.template_defaults["count"] = resX[0]["CountN"].data(); + + if (print_queries) + PrintQuery(query8.str()); + + qlog << query8 << ";\n" << endl; + if (connection.connected()) + { + cout << query8.execute().info() << endl; + ShowWarnings(connection); + Dump(flog, connection, "Spectrum"); + + const auto sec8 = Time().UnixTime()-start8.UnixTime(); + cout << "Execution time: " << sec8 << "s\n"; + + + const mysqlpp::StoreQueryResult res8 = connection.query("SELECT * FROM Spectrum").store(); + + cout << " Bin Flux Error" << endl; + const auto bins = binning_eest.vec(); + for (auto ir=res8.cbegin(); ir!=res8.cend(); ir++) + { + const mysqlpp::Row &row = *ir; + + const uint32_t bin = row[".energy"]; + cout << setw(5) << bins[bin] << ": " << setw(10) << row["Flux"] << " " << setw(10) << row["ErrFlux"] << '\n'; + } + cout << endl; + + // -------------------------------------------------------------------------- + + // Crab Nebula: 1 TeV: 3e-7 / m^2 / s / TeV + // Crab Nebula: 1 TeV: 3e-11 / cm^2 / s / TeV + + sindent indentm(mlog); + + mlog << "void spectrum()\n"; + mlog << "{\n" << indent(4); + mlog << "// Energy Spectrum (e.g. Crab: 3e-11 [cm^-2 s-^1 TeV^-1] @ 1TeV)\n"; + mlog << "TGraphErrors g;\n"; + mlog << "g.SetNameTitle(\"Spectrum\", \"Energy Spectrum\");\n"; + for (auto ir=res8.cbegin(); ir!=res8.cend(); ir++) + { + // This is not efficient but easier to read and safer + const mysqlpp::Row &row = *ir; + + const double hi = row["hi"]; + const double lo = row["lo"]; + const double center = (hi+lo)/2; + const double flux = row["Flux"]; + const double error = row["ErrFlux"]; + + mlog << "g.SetPoint(g.GetN(), " << setw(8) << center << ", " << flux << ");\n"; + mlog << "g.SetPointError(g.GetN()-1, 0, " << error << ");\n"; + } + mlog << "g.GetXaxis()->SetTitle(\"lg(E/TeV)\");\n"; + mlog << "g.GetYaxis()->SetTitle(\"dN/dE [cm^{-2} s^{-1} TeV^{-1}]\");"; + mlog << endl; + } + + WriteHistogram(connection, { + .name = "Signal", + .title = "Signal", + .dir = "Eest/Measurement", + .binningx = binning_eest, + .table = "Spectrum", + .x = ".energy", + .v = "Signal", + .err = "ErrSignal", + .axisx = "Energy lg(E_{est}/GeV)", + .axisy = "Counts" + }); + + WriteHistogram(connection, { + .name = "Background", + .title = "Background", + .dir = "Eest/Measurement", + .binningx = binning_eest, + .table = "Spectrum", + .x = ".energy", + .v = "Background", + .err = "ErrBackground", + .axisx = "Energy lg(E_{est}/GeV)", + .axisy = "Counts" + }); + + WriteHistogram(connection, { + .name = "Excess", + .title = "Excess", + .dir = "Eest/Measurement", + .binningx = binning_eest, + .table = "Spectrum", + .x = ".energy", + .v = "Excess", + .err = "ErrExcess", + .axisx = "Energy lg(E_{est}/GeV)", + .axisy = "Signal - Background (Counts)" + }); + + WriteHistogram(connection, { + .name = "SignalW", + .title = "SignalW", + .dir = "Eest/Simulation/Weighted", + .binningx = binning_eest, + .table = "Spectrum", + .x = ".energy", + .v = "SignalW", + .err = "ErrSignalW", + .axisx = "Energy lg(E_{est}/GeV)", + .axisy = "Weighted" + }); + + WriteHistogram(connection, { + .name = "BackgroundW", + .title = "BackgroundW", + .dir = "Eest/Simulation/Weighted", + .binningx = binning_eest, + .table = "Spectrum", + .x = ".energy", + .v = "BackgroundW", + .err = "ErrBackgroundW", + .axisx = "Energy lg(E_{est}/GeV)", + .axisy = "Weighted" + }); + + WriteHistogram(connection, { + .name = "ExcessW", + .title = "ExcessW", + .dir = "Eest/Simulation/Weighted", + .binningx = binning_eest, + .table = "Spectrum", + .x = ".energy", + .v = "ExcessW", + .err = "ErrExcessW", + .axisx = "Energy lg(E_{est}/GeV)", + .axisy = "Signal - Background (Weighted)" + }); + + WriteHistogram(connection, { + .name = "Significance", + .title = "Significance", + .dir = "Eest/Measurement", + .binningx = binning_eest, + .table = "Spectrum", + .x = ".energy", + .v = "Significance", + .axisx = "Energy lg(E_{est}/GeV)", + .axisy = "Li/Ma Significance" + }); + + WriteHistogram(connection, { + .dir = "Eest", + .name = "Bias", + .title = "Energy Bias", + .binningx = binning_eest, + .table = "Spectrum", + .x = ".energy", + .v = "Bias", + .err = "Resolution", + .axisx = "Energy lg(E_{sim}/GeV)", + .axisy = "lg(E_{est}/E_{sim})", + .stats = false + }); + + WriteHistogram(connection, { + .dir = "Eest", + .name = "Resolution", + .title = "Energy Resolution", + .binningx = binning_eest, + .table = "Spectrum", + .x = ".energy", + .v = "Resolution", + .axisx = "Energy lg(E_{sim}/GeV)", + .axisy = "#sigma(lg(E_{est}/E_{sim}))", + .stats = false + }); + + WriteHistogram(connection, { + .dir = "Eest/Efficiency", + .name = "EfficiencyN", + .title = "Efficiency (Counts)", + .binningx = binning_eest, + .table = "Spectrum", + .x = ".energy", + .v = "EfficiencyN", + .err = "ErrEfficiencyN", + .axisx = "Energy lg(E_{sim}/GeV)", + .axisy = "Efficiency" + }); + + WriteHistogram(connection, { + .dir = "Eest/Efficiency", + .name = "EfficiencyW", + .title = "Efficiency (Weighted)", + .binningx = binning_eest, + .table = "Spectrum", + .x = ".energy", + .v = "EfficiencyW", + .err = "ErrEfficiencyW", + .axisx = "Energy lg(E_{sim}/GeV)", + .axisy = "Efficiency" + }); + + WriteHistogram(connection, { + .name = "Spectrum", + .title = "Differential Energy Spectrum", + .binningx = binning_eest, + .table = "Spectrum", + .x = ".energy", + .v = "Flux", + .err = "ErrFlux", + .axisx = "Energy lg(E/GeV)", + .axisy = "dN/dE [cm^{-2} s^{-1} TeV^{-1}]", + .stats = false + }); + + + // ------------------------------------------------------------------- + // -------------------------- 9: Threshold --------------------------- + + cout << separator("Threshold") << '\n'; + + Time start9; + + /* 02:calculate-threshold.sql */ + + // This query gets the analysis results versus Simulated Energy from the combined table + mysqlpp::Query query9(&connection); + query9 << + // SELECT SUM(CountN) AS CountN, SUM(OnTime) AS OnTime FROM ThetaHist + "CREATE TEMPORARY TABLE Threshold ENGINE=Memory AS\n" + "(\n" + " WITH ThetaSums AS\n" + " (\n" + " SELECT\n" + " SUM(CountN) AS CountSim,\n" + " SUM(OnTime) AS ObsTime\n" + " FROM\n" + " ThetaHist\n" + " ),\n" + " Area AS\n" + " (\n" + " SELECT\n" + " POW(MinImpactHi,2)*PI() AS Area\n" + " FROM\n" + " MonteCarloArea\n" + " ),\n" + " ResultMC AS\n" // Group AnalysisMC by EnergySim Binning + " (\n" + " SELECT\n" + " `.energysim` AS `.energy`,\n" + " ANY_VALUE(ThresholdW) AS ThresholdW,\n" + " ANY_VALUE(ThresholdW2) AS ThresholdW2,\n" + " ANY_VALUE(ThresholdN) AS ThresholdN,\n" + " ANY_VALUE(BiasSim) AS Bias,\n" + " ANY_VALUE(ResolutionSim) AS Resolution\n" + " FROM\n" + " AnalysisMC\n" + " GROUP BY\n" + " `.energy`\n" + " )\n" + " SELECT\n" + " `.energy`, lo, hi,\n" + " ThresholdW,\n" + " SQRT(ThresholdW2) AS ErrThresholdW,\n" + " ThresholdN,\n" + " SQRT(ThresholdN) AS ErrThresholdN,\n" // Scale for Theta-Weights + " ThresholdW * 1000/(POW(10,hi)-POW(10,lo)) / Area / CountSim*ObsTime AS Flux,\n" + " SQRT(ThresholdW2) * 1000/(POW(10,hi)-POW(10,lo)) / Area / CountSim*ObsTime AS ErrFlux,\n" + " Bias,\n" + " Resolution\n" + " FROM\n" + " ResultMC\n" + " INNER JOIN\n" + " BinningEnergySim ON `.energy`=bin\n" + " CROSS JOIN\n" + " ThetaSums, Area\n" + " WHERE\n" + " ThresholdW>0 AND ThresholdW2>0\n" + " ORDER BY\n" + " `.energy`\n" + ")"; + + query9.parse(); + for (auto it=env.cbegin(); it!=env.cend(); it++) + query9.template_defaults[it->first.c_str()] = it->second.c_str(); + + //query9.template_defaults["area"] = area; + //query9.template_defaults["ontime"] = resX[0]["OnTime"].data(); + //query9.template_defaults["count"] = resX[0]["CountN"].data(); + + if (print_queries) + PrintQuery(query9.str()); + + qlog << query9 << ";\n" << endl; + if (connection.connected()) + { + cout << query9.execute().info() << endl; + ShowWarnings(connection); + Dump(flog, connection, "Threshold"); + + const auto sec9 = Time().UnixTime()-start9.UnixTime(); + cout << "Execution time: " << sec9 << "s\n"; + + // -------------------------------------------------------------------------- + + const mysqlpp::StoreQueryResult res9 = connection.query("SELECT * FROM Threshold").store(); + + sindent indentm(mlog, 4); + + mlog << '\n'; + mlog << "// Energy Threshold\n"; + mlog << "TGraphErrors g;\n"; + mlog << "g.SetNameTitle(\"Threshold\", \"Energy Threshold\");\n"; + for (auto ir=res9.cbegin(); ir!=res9.cend(); ir++) + { + // This is not efficient but easier to read and safer + const mysqlpp::Row &row = *ir; + + const double hi = row["hi"]; + const double lo = row["lo"]; + const double center = (hi+lo)/2; + const double width = pow(10, hi)-pow(10, lo); + const double flux = row["Flux"] /width; + const double error = row["ErrFlux"]/width; + + mlog << "g.SetPoint(g.GetN(), " << setw(8) << center << ", " << flux << ");\n"; + mlog << "g.SetPointError(g.GetN()-1, 0, " << error << ");\n"; + } + mlog << "g.GetXaxis()->SetTitle(\"lg(E/TeV)\");\n"; + mlog << "g.GetYaxis()->SetTitle(\"dN/dE [cm^{-2} s^{-1} TeV^{-1}]\");\n"; + mlog << indent(0) << "}" << endl; + } + + WriteHistogram(connection, { + .name = "SimExcessW", + .title = "Weighted Simulated Excess", + .dir = "Esim/Simulation/Weighted", + .binningx = binning_esim, + .table = "Threshold", + .x = ".energy", + .v = "ThresholdW", + .err = "ErrThresholdW", + .axisx = "Energy lg(E_{sim}/GeV)", + .axisy = "Weighted Counts" + }); + + WriteHistogram(connection, { + .name = "SimExcessN", + .title = "Simulated Excess", + .dir = "Esim/Simulation/Counts", + .binningx = binning_esim, + .table = "Threshold", + .x = ".energy", + .v = "ThresholdN", + .err = "ErrThresholdN", + .axisx = "Energy lg(E_{sim}/GeV)", + .axisy = "Counts" + }); + + WriteHistogram(connection, { + .name = "SimSpectrumW", + .title = "Weighted Simulated Excess Spectrum", + .binningx = binning_esim, + .dir = "Esim/Simulation/Weighted", + .table = "Threshold", + .x = ".energy", + .v = "Flux", + .err = "ErrFlux", + .axisx = "Energy lg(E_{sim}/GeV)", + .axisy = "dN/dE [cm^{-2} s^{-1} TeV^{-1}]" + }); + + WriteHistogram(connection, { + .dir = "Esim", + .name = "Bias", + .title = "Energy Bias", + .binningx = binning_esim, + .table = "Threshold", + .x = ".energy", + .v = "Bias", + .err = "Resolution", + .axisx = "Energy lg(E_{sim}/GeV)", + .axisy = "lg(E_{est}/E_{sim})", + .stats = false + }); + + WriteHistogram(connection, { + .dir = "Esim", + .name = "Resolution", + .title = "Energy Resolution", + .binningx = binning_esim, + .table = "Threshold", + .x = ".energy", + .v = "Resolution", + .axisx = "Energy lg(E_{sim}/GeV)", + .axisy = "#sigma(lg(E_{est}/E_{sim}))", + .stats = false + }); + + + // ------------------------------------------------------------------- + // -------------------------- 10: Migration -------------------------- + + cout << separator("Migration") << '\n'; + + Time start10; + + /* 02:obtain-migration-matrix.sql */ + + // This query gets the analysis results versus Estimated Energy from the combined table + mysqlpp::Query query10(&connection); + query10 << + "CREATE TEMPORARY TABLE Migration ENGINE=Memory AS\n" + "(\n" + " SELECT\n" + " `.energyest`,\n" + " `.energysim`,\n" + " BinningEnergySim.lo AS EsimLo,\n" + " BinningEnergySim.hi AS EsimHi,\n" + " BinningEnergyEst.lo AS EestLo,\n" + " BinningEnergyEst.hi AS EestHi,\n" + " ANY_VALUE(MigrationW) AS MigrationW,\n" + " ANY_VALUE(MigrationN) AS MigrationN\n" + // FIXME: Errors + " FROM\n" + " AnalysisMC\n" + " INNER JOIN\n" + " BinningEnergyEst ON `.energyest`=BinningEnergyEst.bin\n" + " INNER JOIN\n" + " BinningEnergySim ON `.energysim`=BinningEnergySim.bin\n" + " GROUP BY\n" + " `.energyest`, `.energysim`\n" + " ORDER BY\n" + " `.energyest`, `.energysim`\n" + ")"; + + if (print_queries) + PrintQuery(query10.str()); + + qlog << query10 << ";\n" << endl; + if (connection.connected()) + { + cout << query10.execute().info() << endl; + ShowWarnings(connection); + Dump(flog, connection, "Migration"); + + const auto sec10 = Time().UnixTime()-start10.UnixTime(); + cout << "Execution time: " << sec10 << "s\n\n"; + } + + WriteHistogram(connection, { + .name = "MigrationN", + .title = "Energy Migration", + .binningx = binning_esim, + .binningy = binning_eest, + .table = "Migration", + .x = ".energysim", + .y = ".energyest", + .v = "MigrationN", + .axisx = "Energy lg(E_{sim}/GeV)", + .axisy = "Energy lg(E_{est}/GeV)", + .axisz = "Counts", + .stats = false + }); + + WriteHistogram(connection, { + .name = "MigrationW", + .title = "Energy Migration", + .binningx = binning_esim, + .binningy = binning_eest, + .table = "Migration", + .x = ".energysim", + .y = ".energyest", + .v = "MigrationW", + .axisx = "Energy lg(E_{sim}/GeV)", + .axisy = "Energy lg(E_{est}/GeV)", + .axisz = "Weighted Counts", + .stats = false + }); + + + + // ------------------------------------------------------------------- + // --------------------------- 11: Summary --------------------------- + + cout << separator("Summary") << '\n'; + const auto sec = Time().UnixTime()-start.UnixTime(); + cout << "Total execution time: " << sec << "s\n" << endl; + + return 0; +}