source: trunk/FACT++/src/spectrum.cc@ 19887

Last change on this file since 19887 was 19887, checked in by tbretz, 6 years ago
In contrast to clang, gcc requires the initilizer list to be in order and complete.
File size: 69.0 KB
Line 
1#include <thread>
2#include <boost/filesystem.hpp>
3#include <boost/range/adaptor/transformed.hpp>
4
5#include "Database.h"
6#include "tools.h"
7#include "Time.h"
8#include "Configuration.h"
9
10#ifdef HAVE_HIGHLIGHT
11#include "srchilite/sourcehighlight.h"
12#include "srchilite/langmap.h"
13#endif
14
15#ifdef HAVE_ROOT
16#include "TFile.h"
17#include "TH1.h"
18#include "TH2.h"
19#endif
20
21using namespace std;
22using boost::adaptors::transformed;
23
24namespace fs = boost::filesystem;
25
26// ------------------------------- Binning ----------------------------------
27
28struct Binning : std::set<double>
29{
30 bool equidist { false };
31
32 Binning() { }
33 Binning(const Binning &m) : std::set<double>(m), equidist(m.equidist) { }
34 Binning(size_t cnt, double mn, double mx) { set(cnt, mn, mx); }
35 void set(size_t cnt, double mn, double mx)
36 {
37 if (cnt==0)
38 return;
39
40 if (empty())
41 equidist = true;
42
43 const double min = ::min(mn, mx);
44 const double max = ::max(mn, mx);
45
46 const double stp = (max-min)/cnt;
47
48 for (size_t i=0; i<=cnt; i++)
49 emplace(min+i*stp);
50 }
51
52 void add(double val)
53 {
54 emplace(val);
55 equidist = false;
56 }
57
58 Binning &operator+=(const vector<double> &v)
59 {
60 if (!v.empty())
61 {
62 insert(v.cbegin(), v.cend());
63 equidist = false;
64 }
65 return *this;
66 }
67
68 Binning operator+(const vector<double> &v)
69 {
70 return Binning(*this) += v;
71 }
72
73 void add(const double &val)
74 {
75 emplace(val);
76 equidist = false;
77 }
78
79 string list() const
80 {
81 return boost::join(*this | transformed([](double d) { return std::to_string(d); }), ",");
82 }
83
84 string str() const
85 {
86 if (!equidist)
87 return list();
88
89 return to_string(size()-1)+":"+to_string(*begin())+","+to_string(*rbegin());
90 }
91
92 vector<double> vec() const { return vector<double>(begin(), end()); }
93
94 //double operator[](size_t i) const { return vec().at(i); }
95};
96
97std::istream &operator>>(std::istream &in, Binning &m)
98{
99 const istreambuf_iterator<char> eos;
100 string txt(istreambuf_iterator<char>(in), eos);
101
102 const boost::regex expr(
103 "([0-9]+)[ /,:;]+"
104 "([-+]?[0-9]*[.]?[0-9]+([eE][-+]?[0-9]+)?)[ /,:;]+"
105 "([-+]?[0-9]*[.]?[0-9]+([eE][-+]?[0-9]+)?)"
106 );
107 boost::smatch match;
108 if (!boost::regex_match(txt, match, expr))
109 throw runtime_error("Could not evaluate binning: "+txt);
110
111 m = Binning(stoi(match[1].str()), stof(match[2].str()), stof(match[4].str()));
112
113 return in;
114}
115
116std::ostream &operator<<(std::ostream &out, const Binning &m)
117{
118 return out << m.str();
119}
120
121// ---------------------------- Configuration -------------------------------
122
123void SetupConfiguration(Configuration &conf)
124{
125 po::options_description control("Calcsource options");
126 control.add_options()
127 ("uri,u", var<string>()
128#if BOOST_VERSION >= 104200
129 ->required()
130#endif
131 , "Database link as in\n\tuser:password@server[:port]/database[?compress=0|1].")
132 ("out,o", var<string>(conf.GetName()), "Defines the prefix (with path) of the output files.")
133 ("dry-run", po_switch(), "Only write the queries to the query.log file. Internally overwrites uri with an empty string.")
134 ;
135
136 po::options_description binnings("Binnings");
137 binnings.add_options()
138 ("theta", var<Binning>(Binning(90, 0, 90)), "Add equidistant bins in theta (degrees). Syntax: N,lo,hi (Number N of equidistant bins between lo and hi)")
139 ("theta-bin", vars<double>(), "Add a bin-edge to the theta binning (degree)")
140 ("esim", var<Binning>(Binning(15, 2, 5)), "Add equidistant bins in log10 simulated energy. Syntax: N,lo,hi (Number N of equidistant bins between lo and hi)")
141 ("esim-bin", vars<double>(), "Add a bin-edge to the binnnig in log10 simulated enegry")
142 ("eest", var<Binning>(Binning(15, 2, 5)), "Add equidistant bins in log10 estimated energy. Syntax: N,lo,hi (Number N of equidistant bins between lo and hi)")
143 ("eest-bin", vars<double>(), "Add a bin-edge to the binning in log10 estimated enegry")
144 ;
145
146 po::options_description queries("Queries");
147 queries.add_options()
148 ("analysis", var<string>("analysis.sql"), "File with the analysis query. A default file is created automatically in the <prefix> directory it does not exist.")
149 ("data", var<string>("data.sql"), "File with the query which creates the data summary. A default file is created automatically in the <prefix> directory it does not exist.")
150 ("simulation", var<string>("simulation.sql"), "File with the query which creates the simulation summary. A default file is created automatically in the <prefix> directory it does not exist.")
151 ;
152
153 po::options_description preparation("Preparation");
154 preparation.add_options()
155 ("source-key", var<uint16_t>(5), "Source key to be used in file selection.")
156 ("selector", vars<string>(), "WHERE clause to be used in file selection.")
157 ("estimator", var<string>()->required(), "Energy estimator to be used.")
158 ("spectrum", var<string>()->required(), "Spectral shape for re-weighting of simulated 'Energy'")
159 ("env.*", var<string>(), "Define a variable that is replaced in all queries automatically.")
160 ;
161
162 po::options_description debug("Debug options");
163 debug.add_options()
164 ("print-connection", po_switch(), "Print database connection information")
165#ifdef HAVE_HIGHLIGHT
166 ("print-queries", po_switch(), "Print all queries to the console. They are automatically highlighted.")
167#else
168 ("print-queries", po_switch(), "Print all queries to the console. (For highlighting recompile with 'libsource-highlight-dev' installed)")
169#endif
170 ("verbose,v", var<uint16_t>(1), "Verbosity (0: quiet, 1: default, 2: more, 3, ...)")
171 ;
172
173 //po::positional_options_description p;
174 //p.add("file", 1); // The 1st positional options (n=1)
175
176 conf.AddOptions(control);
177 conf.AddOptions(binnings);
178 conf.AddOptions(queries);
179 conf.AddOptions(preparation);
180 conf.AddOptions(debug);
181 //conf.SetArgumentPositions(p);
182}
183
184void PrintUsage()
185{
186 //78 .............................................................................
187 cout <<
188 "spectrum - Calculate a spectrum with classical algorithms\n"
189 "\n\n"
190 "Usage: spectrum [-u URI] [options]\n"
191 "\n"
192 ;
193 cout << endl;
194}
195
196
197
198// ----------------------------- Indentation --------------------------------
199
200class sindent : public std::streambuf
201{
202 std::streambuf *sbuf { nullptr };
203 std::ostream *owner { nullptr };
204 int lastch { 0 }; // start-of-line
205 std::string str;
206
207protected:
208 virtual int overflow(int ch)
209 {
210 if (lastch=='\n' && ch != '\n' )
211 sbuf->sputn(str.data(), str.size());
212
213 return sbuf->sputc(lastch = ch);
214 }
215
216public:
217 explicit sindent(std::streambuf *dest, uint16_t indent=0)
218 : sbuf(dest), str(indent, ' ')
219 {
220 }
221
222 explicit sindent(std::ostream& dest, uint16_t indent=0)
223 : sbuf(dest.rdbuf()), owner(&dest), str(indent, ' ')
224 {
225 owner->rdbuf(this);
226 }
227
228 virtual ~sindent()
229 {
230 if (owner)
231 owner->rdbuf(sbuf);
232 }
233
234 void set(uint16_t w) { str.assign(w, ' '); }
235};
236
237struct indent
238{
239 uint16_t w;
240 indent(uint16_t _w=0) : w(_w) { }
241};
242
243std::ostream &operator<<(std::ostream &out, indent m)
244{
245 sindent *buf=dynamic_cast<sindent*>(out.rdbuf());
246 if (buf)
247 buf->set(m.w);
248 return out;
249}
250
251struct separator
252{
253 string s;
254 uint16_t n { 0 };
255 separator(string _s="") : s(_s) { };
256 separator(char c='-', uint16_t _n=78) : s(_n, c), n(_n) { };
257};
258
259std::ostream &operator<<(std::ostream &out, separator m)
260{
261 if (m.n)
262 out << m.s;
263 else
264 {
265 const uint8_t l = (78-m.s.size())/2-3;
266 out << string(l<1?1:l, '-') << "={ " << m.s << " }=" << string(l<1?1:l, '-');
267 if (m.s.size()%2)
268 out << '-';
269 }
270 return out;
271}
272
273// ------------------------------- MySQL++ ----------------------------------
274
275bool ShowWarnings(Database &connection)
276{
277 if (!connection.connected())
278 return true;
279
280 try
281 {
282 const auto resw =
283 connection.query("SHOW WARNINGS").store();
284
285 for (size_t i=0; i<resw.num_rows(); i++)
286 {
287 const mysqlpp::Row &roww = resw[i];
288
289 cout << "\033[31m";
290 cout << roww["Level"] << '[' << roww["Code"] << "]: ";
291 cout << roww["Message"] << "\033[0m" << '\n';
292 }
293 cout << endl;
294 return true;
295
296 }
297 catch (const exception &e)
298 {
299 cerr << "\nSHOW WARNINGS\n\n";
300 cerr << "SQL query failed:\n" << e.what() << '\n' <<endl;
301 return false;
302 }
303}
304
305size_t Dump(ostream &out, Database &connection, const string &table)
306{
307 if (!connection.connected())
308 return 0;
309
310 out << connection.query("SHOW CREATE TABLE `"+table+"`").store()[0][1] << ";\n\n";
311
312 const mysqlpp::StoreQueryResult res = connection.query("SELECT * FROM `"+table+"`").store();
313
314 const string fields = boost::join(res.fields() | transformed([](const mysqlpp::Field &r) { return "`"+string(r.name())+"`"; }), ", ");
315
316 out << "INSERT INTO `" << table << "` ( " << fields << " ) VALUES\n";
317 out << boost::join(res | transformed([](const mysqlpp::Row &row) { return " ( "+boost::join(row | transformed([](const mysqlpp::String &s) { return string(s);}), ", ")+" )";}), ",\n") << ";";
318 out << "\n" << endl;
319
320 return res.num_rows();
321}
322
323void PrintQuery(const string &query)
324{
325#ifdef HAVE_HIGHLIGHT
326 stringstream qstr(query);
327 srchilite::SourceHighlight sourceHighlight("esc256.outlang");
328 sourceHighlight.setStyleFile("esc256.style");
329 sourceHighlight.highlight(qstr, cout, "sql.lang");
330 cout << endl;
331#else
332 cout << query << endl;
333#endif
334}
335
336void CreateBinning(Database &connection, ostream &qlog, const Binning &bins, const string &name)
337{
338 mysqlpp::Query query0(&connection);
339 query0 <<
340 "CREATE TEMPORARY TABLE Binning" << name << "\n"
341 "(\n"
342 " bin SMALLINT UNSIGNED NOT NULL,\n"
343 " lo FLOAT NOT NULL,\n"
344 " hi FLOAT NOT NULL,\n"
345 " PRIMARY KEY (bin) USING HASH\n"
346 ")";
347
348 qlog << query0 << ";\n" << endl;
349 if (connection.connected())
350 query0.execute();
351
352 //connection.query("ALTER TABLE Binning"+name+" AUTO_INCREMENT=0").execute();
353
354 mysqlpp::Query query1(&connection);
355 query1 <<
356 "INSERT INTO Binning" << name << " ( bin, lo, hi ) VALUES\n";
357
358 // FIXME: Case of 1 and 2 bins
359
360 // Bin 0: is the underflow bin...
361 // Bin N: is the overflow bin...
362 // Bin -1: if argument is NULL
363
364 const auto vec = bins.vec();
365 for (size_t i=1; i<vec.size()-2; i++)
366 query1 << " ( " << i << ", " << vec[i-1] << ", " << vec[i] << " ),\n";
367 query1 << " ( " << vec.size()-2 << ", " << vec[vec.size()-2] << ", " << vec[vec.size()-1] << " )\n";
368
369 qlog << query1 << ";\n" << endl;
370
371 if (connection.connected())
372 cout << query1.execute().info() << endl;
373 ShowWarnings(connection);
374}
375
376// ----------------------------- ROOT Histogram -----------------------------
377
378/*
379A bit of hackery, so just sharing for fun.
380
381 #define with(T, ...) ([&]{ T ${}; __VA_ARGS__; return $; }())
382
383And use it like:
384
385 MyFunction(with(Params,
386 $.Name = "Foo Bar",
387 $.Age = 18
388 ));
389
390which expands to:
391
392 MyFunction(([&] {
393 Params ${};
394 $.Name = "Foo Bar", $.Age = 18;
395 return $;
396 }()));
397*/
398struct Histogram
399{
400 // A workaround to be able to set a default also in C++11
401 /*
402 template<typename T, T def>
403 struct Value
404 {
405 T t { def };
406 Value() = default;
407 Value(const T &_t) : t(_t) { }
408 operator T() const { return t; }
409 };*/
410
411 string name;
412 string title;
413 string dir;
414 Binning binningx;
415 Binning binningy;
416 string table;
417 string x;
418 string y;
419 string v;
420 string err;
421 string axisx;
422 string axisy;
423 string axisz;
424 bool stats;
425 //Value<bool,true> stats;
426};
427
428void WriteHistogram(Database &connection, const Histogram &hist)
429{
430#ifdef HAVE_ROOT
431 if (!connection.connected())
432 return;
433
434 mysqlpp::Query query(&connection);
435 query << "SELECT `%0:x` AS X,%1:y `%2:v` AS V%3:err FROM `%4:table`";
436 query.parse();
437
438 query.template_defaults["table"] = hist.table.c_str();
439
440 query.template_defaults["x"] = hist.x.c_str();
441 query.template_defaults["v"] = hist.v.c_str();
442 if (!hist.y.empty())
443 query.template_defaults["y"] = (" `"+hist.y+"` AS Y,").c_str();
444 if (!hist.err.empty())
445 query.template_defaults["err"] = (", `"+hist.err+"` AS E").c_str();
446
447 // PrintQuery(query.str());
448
449 const mysqlpp::StoreQueryResult res = query.store();
450
451 const auto vecx = hist.binningx.vec();
452 const auto vecy = hist.binningy.vec();
453
454 TH1 *h = 0;
455
456 if (hist.y.empty())
457 {
458 h = hist.binningx.equidist ?
459 new TH1D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.front(), vecx.back()) :
460 new TH1D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.data());
461 }
462 else
463 {
464 if (hist.binningy.equidist)
465 {
466 h = hist.binningx.equidist ?
467 new TH2D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.front(), vecx.back(), vecy.size()-1, vecy.front(), vecy.back()) :
468 new TH2D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.data(), vecy.size()-1, vecy.front(), vecy.back());
469 }
470 else
471 {
472 h = hist.binningx.equidist ?
473 new TH2D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.front(), vecx.back(), vecy.size()-1, vecy.front(), vecy.back()) :
474 new TH2D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.data(), vecy.size()-1, vecy.front(), vecy.back());
475 }
476 }
477
478 h->SetXTitle(hist.axisx.c_str());
479 h->SetYTitle(hist.axisy.c_str());
480 h->SetZTitle(hist.axisz.c_str());
481
482 h->SetMarkerStyle(kFullDotMedium);
483
484 h->SetStats(hist.stats);
485
486 for (auto ir=res.cbegin(); ir!=res.cend(); ir++)
487 {
488 const auto &row = *ir;
489
490 const uint32_t x = row["X"];
491 const double v = row["V"];
492 if (x>uint32_t(h->GetNbinsX()+1))
493 continue;
494
495 try
496 {
497 const uint32_t y = row["Y"];
498 if (y>uint32_t(h->GetNbinsY()+1))
499 continue;
500
501 h->SetBinContent(x, y, v);
502
503 }
504 catch (const mysqlpp::BadFieldName &e)
505 {
506 h->SetBinContent(x, v);
507 try
508 {
509 h->SetBinError(x, double(row["E"]));
510 }
511 catch (const mysqlpp::BadFieldName &)
512 {
513 }
514 }
515 }
516
517 gFile->cd();
518 TDirectory *dir = gFile->GetDirectory(hist.dir.c_str());
519 if (dir)
520 dir->cd();
521 else
522 {
523 gFile->mkdir(hist.dir.c_str());
524 gFile->cd(hist.dir.c_str());
525 }
526 h->Write();
527 delete h;
528#endif
529}
530
531// -------------------------------- Resources -------------------------------
532
533#define RESOURCE(TYPE,RC) ([]() { \
534 extern const char _binary_##RC##_start[]; \
535 extern const char _binary_##RC##_end[]; \
536 return TYPE(_binary_##RC##_start, _binary_##RC##_end); \
537})()
538
539string CreateResource(Configuration &conf, const string id, const string def, const string resource)
540{
541 string file = conf.Get<string>(id);
542
543 if (file!=def)
544 {
545 file = conf.GetPrefixedString(id);
546 cout << "Using " << file << "." << endl;
547 return file;
548 }
549
550 file = conf.GetPrefixedString(id);
551
552 cout << "Searching " << file << "... ";
553
554 if (fs::exists(file))
555 {
556 cout << "found." << endl;
557 return file;
558 }
559
560 cout << "creating!" << endl;
561
562 ofstream(file) << resource;
563
564 return file;
565}
566
567// ================================== MAIN ==================================
568
569
570int main(int argc, const char* argv[])
571{
572 Time start;
573
574 Configuration conf(argv[0]);
575 conf.SetPrintUsage(PrintUsage);
576 SetupConfiguration(conf);
577
578 if (!conf.DoParse(argc, argv))
579 return 127;
580
581 // ----------------------------- Evaluate options --------------------------
582
583 const string uri = conf.Get<bool>("dry-run") ? "" : conf.Get<string>("uri");
584 const string out = conf.Get<string>("out");
585 const uint16_t verbose = conf.Get<uint16_t>("verbose");
586
587 const bool print_connection = conf.Get<bool>("print-connection");
588 const bool print_queries = conf.Get<bool>("print-queries");
589
590 Binning binning_theta = conf.Get<Binning>("theta")+conf.Vec<double>("theta-bin");
591 Binning binning_esim = conf.Get<Binning>("esim");
592 Binning binning_eest = conf.Get<Binning>("eest");
593
594 cout << '\n';
595 cout << "Binning 'theta': " << binning_theta.str() << endl;
596 cout << "Binning 'Esim': " << binning_esim.str() << endl;
597 cout << "Binning 'Eest': " << binning_eest.str() << endl;
598
599 const uint16_t source_key = conf.Get<uint16_t>("source-key");
600 const string where = boost::join(conf.Vec<string>("selector"), " AND\n ");
601 const string estimator = conf.Get<string>("estimator");
602 const string spectrum = conf.Get<string>("spectrum");
603 const auto env = conf.GetOptions<string>("env.");
604
605 cout << "\n";
606 const string analysis_sql = CreateResource(conf, "analysis", "analysis.sql", RESOURCE(std::string, spectrum_analysis_sql));
607 const string data_sql = CreateResource(conf, "data", "data.sql", RESOURCE(std::string, spectrum_data_sql));
608 const string simulation_sql = CreateResource(conf, "simulation", "simulation.sql", RESOURCE(std::string, spectrum_simulation_sql));
609 cout << endl;
610
611 const string strzd = binning_theta.list();
612 const string stre = binning_eest.list();
613 const string strth = binning_esim.list();
614
615 // -------------------------------------------------------------------------
616 // Checking for database connection
617
618 Database connection(uri); // Keep alive while fetching rows
619
620 if (!uri.empty())
621 {
622 if (verbose>0)
623 {
624 cout << "Connecting to database...\n";
625 cout << "Client Version: " << mysqlpp::Connection().client_version() << endl;
626 }
627
628 try
629 {
630 connection.connected();
631 }
632 catch (const exception &e)
633 {
634 cerr << "SQL connection failed: " << e.what() << endl;
635 return 1;
636 }
637
638 if (verbose>0)
639 {
640 cout << "Server Version: " << connection.server_version() << endl;
641 cout << "Connected to " << connection.ipc_info() << endl;
642 }
643
644 if (print_connection)
645 {
646 try
647 {
648 const auto res1 = connection.query("SHOW STATUS LIKE 'Compression'").store();
649 cout << "Compression of database connection is " << string(res1[0][1]) << endl;
650
651 const auto res2 = connection.query("SHOW STATUS LIKE 'Ssl_cipher'").store();
652 cout << "Connection to databases is " << (string(res2[0][1]).empty()?"UNENCRYPTED":"ENCRYPTED ("+string(res2[0][1])+")") << endl;
653 }
654 catch (const exception &e)
655 {
656 cerr << "\nSHOW STATUS LIKE 'Compression'\n\n";
657 cerr << "SQL query failed:\n" << e.what() << endl;
658 return 9;
659 }
660 }
661 }
662
663 // -------------------------------------------------------------------------
664 // Create log streams
665
666 ofstream qlog(out+".query.sql");
667 ofstream flog(connection.connected() ? out+".dump.sql" : "");
668 ofstream mlog(connection.connected() ? out+".C" : "");
669
670 cout << "\n";
671 cout << "Queries will be logged to " << out << ".query.sql\n";
672 if (connection.connected())
673 {
674 cout << "Tables will be dumped to " << out << ".dump.sql\n";
675 cout << "ROOT macro will be written to " << out << ".C\n";
676 }
677
678#ifdef HAVE_ROOT
679 TFile root(connection.connected() ? (out+".hist.root").c_str() : "", "RECREATE");
680 if (connection.connected())
681 {
682 if (root.IsZombie())
683 return 10;
684 cout << "Histograms will be written to " << out << ".hist.root\n";
685 }
686#endif
687
688 cout << endl;
689
690 // FIMXE: Implement SYNTAX check on spectrum, estimator and selector
691
692 // -------------------------------------------------------------------
693 // --------------------------- 0: Binnings ---------------------------
694
695 cout << separator("Binnings") << '\n';
696
697 CreateBinning(connection, qlog, binning_theta, "Theta");
698 CreateBinning(connection, qlog, binning_eest, "EnergyEst");
699 CreateBinning(connection, qlog, binning_esim, "EnergySim");
700
701 Dump(flog, connection, "BinningTheta");
702 Dump(flog, connection, "BinningEnergyEst");
703 Dump(flog, connection, "BinningEnergySim");
704
705 // -------------------------------------------------------------------
706 // --------------------------- 1: DataFiles --------------------------
707 cout << separator("DataFiles") << '\n';
708
709 Time start1;
710
711 /* 01:get-data-files.sql */
712 mysqlpp::Query query1(&connection);
713 query1 <<
714 "CREATE TEMPORARY TABLE DataFiles\n"
715 "(\n"
716 " FileId INT UNSIGNED NOT NULL,\n"
717 " PRIMARY KEY (FileId) USING HASH\n"
718 ") ENGINE=Memory\n"
719 "AS\n"
720 "(\n"
721 " SELECT\n"
722 " FileId\n"
723 " FROM\n"
724 " factdata.RunInfo\n"
725 " WHERE\n"
726 " fRunTypeKEY=1 AND fSourceKEY=%100:source AND\n"
727 " %101:where\n"
728 " ORDER BY\n"
729 " FileId\n" // In order: faster
730 ")";
731
732 query1.parse();
733 for (auto it=env.cbegin(); it!=env.cend(); it++)
734 query1.template_defaults[it->first.c_str()] = it->second.c_str();
735
736 query1.template_defaults["source"] = to_string(source_key).c_str();
737 query1.template_defaults["where"] = where.c_str();
738
739 if (print_queries)
740 PrintQuery(query1.str());
741
742 qlog << query1 << ";\n" << endl;
743 if (connection.connected())
744 {
745 cout << query1.execute().info() << endl;
746 ShowWarnings(connection);
747 Dump(flog, connection, "DataFiles");
748
749 const auto sec1 = Time().UnixTime()-start1.UnixTime();
750 cout << "Execution time: " << sec1 << "s\n" << endl;
751 }
752
753 // FIXME: Setup Zd binning depending on Data
754
755 // -------------------------------------------------------------------
756 // ------------------------ 2: ObservationTime -----------------------
757
758 cout << separator("ObservationTime") << '\n';
759
760 Time start2;
761
762 /* 02:get-observation-time.sql */
763 mysqlpp::Query query2(&connection);
764 query2 <<
765 "CREATE TEMPORARY TABLE ObservationTime\n"
766 "(\n"
767 " `.theta` SMALLINT UNSIGNED NOT NULL,\n"
768 " OnTime FLOAT NOT NULL,\n"
769 " PRIMARY KEY (`.theta`) USING HASH\n"
770 ") ENGINE=Memory\n"
771 "AS\n"
772 "(\n"
773 " SELECT\n"
774 " INTERVAL(fZenithDistanceMean, %100:bins) AS `.theta`,\n"
775 " SUM(TIME_TO_SEC(TIMEDIFF(fRunStop,fRunStart))*fEffectiveOn) AS OnTime\n"
776 " FROM\n"
777 " DataFiles\n"
778 " LEFT JOIN \n"
779 " factdata.RunInfo USING (FileId)\n"
780 " GROUP BY\n"
781 " `.theta`\n"
782 " ORDER BY\n"
783 " `.theta`\n"
784 ")";
785
786 query2.parse();
787 for (auto it=env.cbegin(); it!=env.cend(); it++)
788 query2.template_defaults[it->first.c_str()] = it->second.c_str();
789
790 query2.template_defaults["bins"] = strzd.c_str();
791
792 if (print_queries)
793 PrintQuery(query2.str());
794
795 qlog << query2 << ";\n" << endl;
796 if (connection.connected())
797 {
798 cout << query2.execute().info() << endl;
799 ShowWarnings(connection);
800 Dump(flog, connection, "ObservationTime");
801
802 const auto sec2 = Time().UnixTime()-start2.UnixTime();
803 cout << "Execution time: " << sec2 << "s\n\n";
804 }
805
806 /*
807
808 SELECT
809 min(first_number) as first_number,
810 max(last_number) as last_number
811 from
812 (
813 select
814 first_number,
815 last_number
816 sum(grp) over(order by first_number) as grp
817 from
818 (
819 select
820 bin AS first_number,
821 LEAD(bin) AS last_number
822 IF(first_number <> lag(bin, 2) over (order by bin) + 1, 1, 0) as grp
823 from
824 t1
825 )
826 )
827 group by grp
828 order by 1
829 */
830
831 /*
832WITH Hallo AS (WITH MyTest AS (WITH MyTable AS (SELECT bin FROM Test
833WHERE cnt>0 GROUP BY bin ORDER BY bin) SELECT bin AS first_number,
834LAG(bin) OVER(ORDER BY bin) AS prev, LEAD(bin) OVER(ORDER BY bin) AS
835last_number, IF(bin <> (LAG(bin) OVER(ORDER BY bin)+1), 1, 0) AS
836grp_prev, IF(bin <> (LEAD(bin) OVER(ORDER BY bin)-1), 1, 0) as grp_nxt
837FROM MyTable ) SELECT first_number, last_number, sum(grp_nxt)
838over(order by first_number) as grp FROM MyTest) SELECT * FROM Hallo
839
840
841 WITH MyTable AS
842(SELECT bin FROM Test WHERE cnt>0 GROUP BY bin ORDER BY bin)
843SELECT bin, LAG(bin) OVER(ORDER BY bin) AS prev, LEAD(bin) OVER(ORDER BY bin) AS nxt,
844IF(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
845
846"-1","0"
847"0","0"
848"2","3" 2-4
849"3","5"
850"4","6"
851"5","0"
852"6","0"
853"7","7" 7-7
854"8","0"
855"9","9" 9-12
856"10","10"
857"11","11"
858"12","12"
859"13","0"
860"14","0"
861"15","0"
862"16","16" 16-19
863"17","17"
864"18","18"
865"19","19"
866 */
867
868 // -------------------------------------------------------------------
869 // ------------------------ 3: MonteCarloFiles -----------------------
870
871 cout << separator("MonteCarloFiles") << '\n';
872
873 Time start3;
874
875 /* 02:get-monte-carlo-files.sql */
876 mysqlpp::Query query3(&connection);
877 query3 <<
878 "CREATE TEMPORARY TABLE MonteCarloFiles\n"
879 "(\n"
880 " FileId INT UNSIGNED NOT NULL,\n"
881 " PRIMARY KEY (FileId) USING HASH\n"
882 ") ENGINE=Memory\n"
883 "AS\n"
884 "(\n"
885 " SELECT\n"
886 " FileId\n"
887 " FROM\n"
888 " ObservationTime\n"
889 " LEFT JOIN\n"
890 " BinningTheta ON `.theta`=bin\n"
891 " LEFT JOIN\n"
892 " factmc.RunInfoMC\n"
893 " ON\n"
894 //" ThetaMin BETWEEN lo AND hi OR ThetaMax BETWEEN lo AND hi\n" // Includes BOTH edges
895 " (ThetaMin>=lo AND ThetaMin<hi) OR (ThetaMax>lo AND ThetaMax<=hi)\n"
896 " WHERE\n"
897 " PartId=1 AND\n"
898 " FileId%%2=0\n"
899 " ORDER BY\n"
900 " FileId\n" // In order: faster
901 ")";
902/*
903 " SELECT\n"
904 " FileId\n"
905 " FROM\n"
906 " factmc.RunInfoMC\n"
907 " WHERE\n"
908 " PartId=1 AND\n"
909 " FileId%%2=0 AND\n"
910 " (%0:range)\n"
911 ")";
912*/
913
914 query3.parse();
915 for (auto it=env.cbegin(); it!=env.cend(); it++)
916 query3.template_defaults[it->first.c_str()] = it->second.c_str();
917
918 if (print_queries)
919 PrintQuery(query3.str());
920
921 qlog << query3 << ";\n" << endl;
922 if (connection.connected())
923 {
924 cout << query3.execute().info() << endl;
925 ShowWarnings(connection);
926 Dump(flog, connection, "MonteCarloFiles");
927
928 const auto sec3 = Time().UnixTime()-start3.UnixTime();
929 cout << "Execution time: " << sec3 << "s\n\n";
930 }
931
932 // -------------------------------------------------------------------
933 // ----------------------- 3b: Monte Carlo Area ----------------------
934
935 cout << separator("MonteCarloArea") << '\n';
936
937 Time start3b;
938
939 /* 02:get-monte-carlo-files.sql */
940 mysqlpp::Query query3b(&connection);
941 query3b <<
942 "CREATE TEMPORARY TABLE MonteCarloArea ENGINE=Memory\n"
943 "AS\n"
944 "(\n"
945 " SELECT\n"
946 " MIN(`CSCAT[1]`) AS MinImpactLo,\n"
947 " MAX(`CSCAT[1]`) AS MaxImpactLo,\n"
948 " MIN(`CSCAT[2]`) AS MinImpactHi,\n"
949 " MAX(`CSCAT[2]`) AS MaxImpactHi\n"
950 " FROM\n"
951 " MonteCarloFiles\n"
952 " LEFT JOIN\n"
953 " factmc.CorsikaSetup ON FileId=RUNNR\n"
954 " GROUP BY\n"
955 " `CSCAT[1]`, `CSCAT[2]`\n"
956 " ORDER BY\n"
957 " MaxImpactHi\n"
958 ")";
959
960 query3b.parse();
961 for (auto it=env.cbegin(); it!=env.cend(); it++)
962 query3b.template_defaults[it->first.c_str()] = it->second.c_str();
963
964 if (print_queries)
965 PrintQuery(query3b.str());
966
967 qlog << query3b << ";\n" << endl;
968 if (connection.connected())
969 {
970 cout << query3b.execute().info() << endl;
971 ShowWarnings(connection);
972 if (Dump(flog, connection, "MonteCarloArea")!=1)
973 {
974 cerr << "Impact range inconsistent!" << endl;
975 return 1;
976 }
977
978 const auto sec3b = Time().UnixTime()-start3b.UnixTime();
979 cout << "Execution time: " << sec3b << "s\n\n";
980 }
981
982
983 // -------------------------------------------------------------------
984 // -------------------------- 4: ThetaHist ---------------------------
985
986 cout << separator("ThetaHist") << '\n';
987
988 Time start4;
989
990 /* 02:get-theta-distribution.sql */
991 mysqlpp::Query query4(&connection);
992 query4 <<
993 "CREATE TEMPORARY TABLE ThetaHist\n"
994 "(\n"
995 " `.theta` SMALLINT UNSIGNED NOT NULL,\n"
996 " lo DOUBLE NOT NULL COMMENT 'Lower edge of zenith distance bin in degree',\n"
997 " hi DOUBLE NOT NULL COMMENT 'Upper edge of zenith distance bin in degree',\n"
998 " CountN INT UNSIGNED NOT NULL,\n"
999 " ErrCountN DOUBLE NOT NULL,\n"
1000 " OnTime FLOAT NOT NULL,\n"
1001 " ZdWeight DOUBLE NOT NULL COMMENT 'tau(delta theta)',\n"
1002 " ErrZdWeight DOUBLE NOT NULL COMMENT 'sigma(tau)',\n"
1003 " PRIMARY KEY (`.theta`) USING HASH\n"
1004 ") ENGINE=Memory\n"
1005 "AS\n"
1006 "(\n"
1007 " WITH EventCount AS\n"
1008 " (\n"
1009 " SELECT\n"
1010 " INTERVAL(DEGREES(Theta), %100:bins) AS `.theta`,\n"
1011 " COUNT(*) AS CountN\n"
1012 " FROM\n"
1013 " MonteCarloFiles\n"
1014 " LEFT JOIN\n"
1015 " factmc.OriginalMC USING(FileId)\n"
1016 " GROUP BY\n"
1017 " `.theta`\n"
1018 " )\n"
1019 " SELECT\n"
1020 " `.theta`, lo, hi,\n"
1021 " CountN,\n"
1022 " SQRT(CountN) AS ErrCountN,\n"
1023 " OnTime,\n"
1024 " OnTime/CountN AS ZdWeight,\n"
1025 " (OnTime/CountN)*SQRT(POW(1/300, 2) + 1/CountN) AS ErrZdWeight\n"
1026 " FROM\n"
1027 " ObservationTime\n"
1028 " LEFT JOIN\n"
1029 " EventCount USING(`.theta`)\n"
1030 " LEFT JOIN\n"
1031 " BinningTheta ON `.theta`=bin\n"
1032 " ORDER BY\n"
1033 " `.theta`\n"
1034 ")";
1035
1036 query4.parse();
1037 for (auto it=env.cbegin(); it!=env.cend(); it++)
1038 query4.template_defaults[it->first.c_str()] = it->second.c_str();
1039
1040 query4.template_defaults["bins"] = strzd.c_str();
1041
1042 if (print_queries)
1043 PrintQuery(query4.str());
1044
1045 qlog << query4 << ";\n" << endl;
1046 if (connection.connected())
1047 {
1048 cout << query4.execute().info() << endl;
1049 ShowWarnings(connection);
1050 Dump(flog, connection, "ThetaHist");
1051
1052 const auto sec4 = Time().UnixTime()-start4.UnixTime();
1053 cout << "Execution time: " << sec4 << "s\n";
1054
1055
1056 if (verbose>0)
1057 {
1058 const mysqlpp::StoreQueryResult res4 = connection.query("SELECT * FROM ThetaHist").store();
1059
1060 cout << " Bin MC Counts OnTime ZdWeight\n";
1061 const auto bins = binning_theta.vec();
1062 for (auto ir=res4.cbegin(); ir!=res4.cend(); ir++)
1063 {
1064 const mysqlpp::Row &row = *ir;
1065
1066 const uint32_t bin = row[".theta"];
1067 cout << setw(5) << bins[bin] << ": " << setw(10) << row["CountN"] << " " << setw(10) << row["OnTime"] << " " << setw(10) << row["ZdWeight"] << '\n';
1068 }
1069 cout << endl;
1070 }
1071 }
1072
1073 WriteHistogram(connection, {
1074 .name = "OnTime",
1075 .title = "Effective on time",
1076 .dir = "Zd",
1077 .binningx = binning_theta,
1078 .binningy = {},
1079 .table = "ThetaHist",
1080 .x = ".theta",
1081 .y = "",
1082 .v = "OnTime",
1083 .err = "",
1084 .axisx = "Zenith Distance #theta [#circ]",
1085 .axisy = "Eff. on time [s]",
1086 .axisz = "",
1087 .stats = true
1088 });
1089
1090 WriteHistogram(connection, {
1091 .name = "CountN",
1092 .title = "Simulated Zenith Distance",
1093 .dir = "Zd",
1094 .binningx = binning_theta,
1095 .binningy = {},
1096 .table = "ThetaHist",
1097 .x = ".theta",
1098 .y = "",
1099 .v = "CountN",
1100 .err = "ErrCountN",
1101 .axisx = "Zenith Distance #theta [#circ]",
1102 .axisy = "Counts",
1103 .axisz = "",
1104 .stats = true
1105 });
1106
1107 WriteHistogram(connection, {
1108 .name = "ZdWeight",
1109 .title = "Zenith Distance Weight",
1110 .dir = "Zd",
1111 .binningx = binning_theta,
1112 .binningy = {},
1113 .table = "ThetaHist",
1114 .x = ".theta",
1115 .y = "",
1116 .v = "ZdWeight",
1117 .err = "ErrZdWeight",
1118 .axisx = "Zenith Distance #theta [#circ]",
1119 .axisy = "Weight [s]",
1120 .axisz = "",
1121 .stats = true
1122 });
1123
1124
1125 // -------------------------------------------------------------------
1126 // ------------------------- 5: AnalysisData -------------------------
1127
1128 cout << separator("AnalysisData") << '\n';
1129
1130 Time start5;
1131
1132 /* 02:analyze-data.sql */
1133 mysqlpp::Query query5(&connection);
1134 sindent indent5(query5);
1135 query5 <<
1136 "CREATE TEMPORARY TABLE AnalysisData\n"
1137 "(\n"
1138 " `.energy` SMALLINT UNSIGNED NOT NULL,\n"
1139 " `Signal` DOUBLE NOT NULL,\n"
1140 " `Background` DOUBLE NOT NULL,\n"
1141 " `Excess` DOUBLE NOT NULL,\n"
1142 " `Significance` DOUBLE NOT NULL,\n"
1143 " `ErrExcess` DOUBLE NOT NULL,\n"
1144 " PRIMARY KEY (`.energy`) USING HASH\n"
1145 ") ENGINE=Memory\n"
1146 "AS\n"
1147 "(\n"
1148 " WITH Excess AS\n"
1149 " (\n" << indent(6)
1150 << ifstream(analysis_sql).rdbuf() << indent(0) <<
1151 " ),\n" <<
1152 " Result AS\n"
1153 " (\n" << indent(6)
1154 << ifstream(data_sql).rdbuf() << indent(0) << // Must end with EOL and not in the middle of a comment
1155 " )\n"
1156 " SELECT * FROM Result\n"
1157 ")";
1158
1159 query5.parse();
1160 for (auto it=env.cbegin(); it!=env.cend(); it++)
1161 query5.template_defaults[it->first.c_str()] = it->second.c_str();
1162
1163 //query5.template_defaults["columns"] = "FileId, EvtNumber,";
1164 query5.template_defaults["columns"] = "";
1165 query5.template_defaults["files"] = "DataFiles";
1166 query5.template_defaults["runinfo"] = "factdata.RunInfo";
1167 query5.template_defaults["events"] = "factdata.Images";
1168 query5.template_defaults["positions"] = "factdata.Position";
1169
1170 query5.template_defaults["bins"] = stre.c_str();
1171 query5.template_defaults["estimator"] = estimator.c_str();
1172
1173 if (print_queries)
1174 PrintQuery(query5.str());
1175
1176 qlog << query5 << ";\n" << endl;
1177 if (connection.connected())
1178 {
1179 cout << query5.execute().info() << endl;
1180 ShowWarnings(connection);
1181 Dump(flog, connection, "AnalysisData");
1182
1183 const auto sec5 = Time().UnixTime()-start5.UnixTime();
1184 cout << "Execution time: " << sec5 << "s\n";
1185
1186 if (verbose>0)
1187 {
1188 const mysqlpp::StoreQueryResult res5 = connection.query("SELECT * FROM AnalysisData").store();
1189
1190 cout << " Bin Signal Background Excess Significance Error" << endl;
1191 for (auto row=res5.cbegin(); row!=res5.cend(); row++)
1192 {
1193 for (auto it=row->begin(); it!=row->end(); it++)
1194 cout << setw(10) << *it << " ";
1195 cout << '\n';
1196 }
1197 cout << endl;
1198 }
1199 }
1200
1201
1202 // -------------------------------------------------------------------
1203 // --------------------------- 6: ResultMC ---------------------------
1204
1205 cout << separator("ResultMC") << '\n';
1206
1207 Time start6;
1208
1209 /* 02:analyze-simulation.sql */
1210
1211 // This table combines the analysis results vs. Binning in Estimated Energy and Simulated Energy
1212 mysqlpp::Query query6(&connection);
1213 sindent indent6(query6);
1214 query6 <<
1215 "CREATE TEMPORARY TABLE AnalysisMC\n"
1216 "(\n"
1217 " `.energyest` SMALLINT UNSIGNED NOT NULL,\n"
1218 " `.energysim` SMALLINT UNSIGNED NOT NULL,\n"
1219 " SignalW DOUBLE NOT NULL,\n"
1220 " BackgroundW DOUBLE NOT NULL,\n"
1221 " ThresholdW DOUBLE NOT NULL,\n"
1222 " SignalN DOUBLE NOT NULL,\n"
1223 " BackgroundN DOUBLE NOT NULL,\n"
1224 " ThresholdN DOUBLE NOT NULL,\n"
1225 " ExcessW DOUBLE NOT NULL,\n"
1226 " ExcessN DOUBLE NOT NULL,\n"
1227 " ErrExcessW DOUBLE NOT NULL,\n"
1228 " ErrExcessN DOUBLE NOT NULL,\n"
1229 " BiasEst DOUBLE NOT NULL,\n"
1230 " BiasSim DOUBLE NOT NULL,\n"
1231 " ResolutionEst DOUBLE,\n"
1232 " ResolutionSim DOUBLE,\n"
1233 " INDEX (`.energyest`) USING HASH,\n"
1234 " INDEX (`.energysim`) USING HASH\n"
1235 ") ENGINE=Memory\n"
1236 "AS\n"
1237 "(\n"
1238 " WITH Excess AS\n"
1239 " (\n" << indent(6)
1240 << ifstream(analysis_sql).rdbuf() << indent(0) <<
1241 " ),\n" <<
1242 " Result AS\n"
1243 " (\n" << indent(6)
1244 << ifstream(simulation_sql).rdbuf() << indent(0) << // Must end with EOL and not in the middle of a comment
1245 " )\n"
1246 " SELECT * FROM Result\n"
1247 ")";
1248
1249 query6.parse();
1250 for (auto it=env.cbegin(); it!=env.cend(); it++)
1251 query6.template_defaults[it->first.c_str()] = it->second.c_str();
1252
1253 //query6.template_defaults["columns"] = "FileId, EvtNumber, CorsikaNumReuse,";
1254 query6.template_defaults["columns"] = "Zd, Energy, SpectralIndex,";
1255 query6.template_defaults["files"] = "MonteCarloFiles";
1256 query6.template_defaults["runinfo"] = "factmc.RunInfoMC";
1257 query6.template_defaults["events"] = "factmc.EventsMC";
1258 query6.template_defaults["positions"] = "factmc.PositionMC";
1259
1260 query6.template_defaults["energyest"] = stre.c_str();
1261 query6.template_defaults["energysim"] = strth.c_str();
1262 query6.template_defaults["theta"] = strzd.c_str();
1263 query6.template_defaults["spectrum"] = spectrum.c_str();
1264 query6.template_defaults["estimator"] = estimator.c_str();
1265
1266 if (print_queries)
1267 PrintQuery(query6.str());
1268
1269 qlog << query6 << ";\n" << endl;
1270 if (connection.connected())
1271 {
1272 cout << query6.execute().info() << endl;
1273 ShowWarnings(connection);
1274 Dump(flog, connection, "AnalysisMC");
1275
1276 const auto sec6 = Time().UnixTime()-start6.UnixTime();
1277 cout << "Execution time: " << sec6 << "s\n\n";
1278 }
1279
1280
1281 // -------------------------------------------------------------------
1282 // ----------------------- 7: SimulatedSpectrum ----------------------
1283
1284 cout << separator("SimulatedSpectrum") << '\n';
1285
1286 // [lo^(1+gammaa) - hi^(1+gamma)] / (1+gamma)
1287
1288 Time start7;
1289 /* 02:get-corsika-events.sql */
1290
1291 // FIXME: Theta weights?
1292 // FIXME: energysim binning
1293 mysqlpp::Query query7(&connection);
1294 query7 <<
1295 "CREATE TEMPORARY TABLE SimulatedSpectrum\n"
1296 "(\n"
1297 " `.energy` SMALLINT UNSIGNED NOT NULL COMMENT 'Bin Index [MC Energy]',\n"
1298 " CountN DOUBLE NOT NULL,\n" // COMMENT 'Number of events',\n"
1299 " CountW DOUBLE NOT NULL,\n" // COMMENT 'Sum of weights, reweighted sum',\n"
1300 " CountW2 DOUBLE NOT NULL,\n" // COMMENT 'Sum of square of weights'\n"
1301 " PRIMARY KEY (`.energy`) USING HASH\n"
1302 ") ENGINE=Memory\n"
1303 "AS\n"
1304 "(\n"
1305 " SELECT\n"
1306 " INTERVAL(LOG10(Energy), %100:energyest) AS `.energy`,\n"
1307 " COUNT(*) AS CountN,\n"
1308 " SUM( (%101:spectrum)/pow(Energy, SpectralIndex) ) AS CountW,\n" // Contents is: CountW
1309 " SUM(POW((%101:spectrum)/pow(Energy, SpectralIndex),2)) AS CountW2\n" // Error is: SQRT(CountW2)
1310 " FROM\n"
1311 " MonteCarloFiles\n"
1312 " LEFT JOIN\n"
1313 " factmc.RunInfoMC USING (FileId)\n"
1314 " LEFT JOIN\n"
1315 " factmc.OriginalMC USING (FileId)\n"
1316 " GROUP BY\n"
1317 " `.energy`\n"
1318 " ORDER BY\n"
1319 " `.energy`\n"
1320 ")";
1321
1322 query7.parse();
1323 for (auto it=env.cbegin(); it!=env.cend(); it++)
1324 query7.template_defaults[it->first.c_str()] = it->second.c_str();
1325
1326 //query7.template_defaults["list"] = listmc.c_str();
1327 query7.template_defaults["energyest"] = stre.c_str();
1328 query7.template_defaults["spectrum"] = spectrum.c_str();
1329
1330 if (print_queries)
1331 PrintQuery(query7.str());
1332
1333 qlog << query7 << ";\n" << endl;
1334 if (connection.connected())
1335 {
1336 cout << query7.execute().info() << endl;
1337 ShowWarnings(connection);
1338 Dump(flog, connection, "SimulatedSpectrum");
1339
1340 const auto sec7 = Time().UnixTime()-start7.UnixTime();
1341 cout << "Execution time: " << sec7 << "s\n";
1342
1343
1344 if (verbose>0)
1345 {
1346 const mysqlpp::StoreQueryResult res7 = connection.query("SELECT * FROM SimulatedSpectrum").store();
1347
1348 cout << " Bin CountW CountW2" << endl;
1349 const auto bins = binning_esim.vec();
1350 for (auto ir=res7.cbegin(); ir!=res7.cend(); ir++)
1351 {
1352 const mysqlpp::Row &row = *ir;
1353
1354 const uint32_t bin = row[".energy"];
1355 cout << setw(5) << bins[bin] << ": " << setw(10) << row["CountW"] << " " << setw(10) << row["CountW2"] << '\n';
1356 }
1357 cout << endl;
1358 }
1359 }
1360
1361
1362 // -------------------------------------------------------------------
1363 // --------------------------- 8: Spectrum ---------------------------
1364
1365 cout << separator("Spectrum") << '\n';
1366
1367 Time start8;
1368
1369 /* 02:calculate-spectrum.sql */
1370
1371 mysqlpp::Query query8(&connection);
1372 query8 <<
1373 // SELECT SUM(CountN) AS CountN, SUM(OnTime) AS OnTime FROM ThetaHist
1374 "CREATE TEMPORARY TABLE Spectrum\n"
1375 "(\n"
1376 " `.energy` SMALLINT UNSIGNED NOT NULL COMMENT 'Bin Index [Energy]' PRIMARY KEY,\n"
1377 " lo DOUBLE NOT NULL COMMENT 'Lower edge of energy bin in lg(E/GeV)',\n"
1378 " hi DOUBLE NOT NULL COMMENT 'Upper edge of energy bin in lg(E/GeV)',\n"
1379 " `Signal` DOUBLE NOT NULL COMMENT 'Number of signal events',\n"
1380 " `Background` DOUBLE NOT NULL COMMENT 'Average number of background events',\n"
1381 " `Excess` DOUBLE NOT NULL COMMENT 'Number of excess events',\n"
1382 " ErrSignal DOUBLE NOT NULL COMMENT 'Poisson error on number of signal events',\n"
1383 " ErrBackground DOUBLE NOT NULL COMMENT 'Poisson error on number of background events',\n"
1384 " `ErrExcess` DOUBLE NOT NULL COMMENT 'Error of excess events',\n"
1385 " `Significance` DOUBLE NOT NULL COMMENT 'Li/Ma sigficance',\n"
1386 " `ExcessN` DOUBLE NOT NULL COMMENT 'Number of excess events in simulated data',\n"
1387 " `ExcessW` DOUBLE NOT NULL COMMENT 'Weighted number of excess events in simulated data',\n"
1388 " `ErrExcessN` DOUBLE NOT NULL COMMENT 'Error or number of excess events in simulated data',\n"
1389 " `ErrExcessW` DOUBLE NOT NULL COMMENT 'Error of weighted number of excess events in simulated data',\n"
1390 " SignalW DOUBLE NOT NULL COMMENT 'Weighted number of signal events in simulated data',\n"
1391 " BackgroundW DOUBLE NOT NULL COMMENT 'Weighted number of background events in simulated data',\n"
1392 " ErrSignalW DOUBLE NOT NULL COMMENT 'Error of weighted number of signal events in simulated data',\n"
1393 " ErrBackgroundW DOUBLE NOT NULL COMMENT 'Error of weighted number of background events in simulated data',\n"
1394 " Flux DOUBLE NOT NULL COMMENT 'dN/dA/dt [cm^-2 s-^1 TeV^-1]',\n"
1395 " ErrFlux DOUBLE NOT NULL COMMENT 'dN/dA/dt [cm^-2 s-^1 TeV^-1]',\n"
1396 " Bias DOUBLE NOT NULL COMMENT 'Energy Bias, average residual in lg(E)',\n"
1397 " Resolution DOUBLE NOT NULL COMMENT 'Energy resolution, standard divation of residual in lg(E)',\n"
1398 " EfficiencyN DOUBLE NOT NULL COMMENT 'Simulated cut efficiency (weighted)',\n"
1399 " EfficiencyW DOUBLE NOT NULL COMMENT 'Simulated cut efficiency (unweighted)',\n"
1400 " ErrEfficiencyN DOUBLE NOT NULL COMMENT 'Error of simulated cut efficiency (weighted)',\n"
1401 " ErrEfficiencyW DOUBLE NOT NULL COMMENT 'Error of simulated cut efficiency (unweighted)'\n"
1402 ") ENGINE=Memory\n"
1403 "AS\n"
1404 "(\n"
1405 " WITH ThetaSums AS\n"
1406 " (\n"
1407 " SELECT\n"
1408 " SUM(CountN) AS CountSim,\n"
1409 " SUM(OnTime) AS ObsTime\n"
1410 " FROM\n"
1411 " ThetaHist\n"
1412 " ),\n"
1413 " Area AS\n"
1414 " (\n"
1415 " SELECT\n"
1416 " POW(MinImpactHi,2)*PI() AS Area\n"
1417 " FROM\n"
1418 " MonteCarloArea\n"
1419 " ),\n"
1420 " ResultMC AS\n" // Group AnalysisMC by EnergyEst Binning
1421 " (\n"
1422 " SELECT\n"
1423 " `.energyest` AS `.energy`,\n"
1424 " ANY_VALUE(SignalW) AS SignalW,\n"
1425 " ANY_VALUE(SignalW2) AS SignalW2,\n"
1426 " ANY_VALUE(BackgroundW) AS BackgroundW,\n"
1427 " ANY_VALUE(BackgroundW2) AS BackgroundW2,\n"
1428 " ANY_VALUE(SignalN) AS SignalN,\n"
1429 " ANY_VALUE(BackgroundN) AS BackgroundN,\n"
1430 " ANY_VALUE(ExcessW) AS ExcessW,\n"
1431 " ANY_VALUE(ExcessN) AS ExcessN,\n"
1432 " ANY_VALUE(ErrExcessW) AS ErrExcessW,\n"
1433 " ANY_VALUE(ErrExcessN) AS ErrExcessN,\n"
1434 " ANY_VALUE(BiasEst) AS Bias,\n"
1435 " ANY_VALUE(ResolutionEst) AS Resolution\n"
1436 " FROM\n"
1437 " AnalysisMC\n"
1438 " GROUP BY\n"
1439 " `.energy`\n"
1440 " ORDER BY\n"
1441 " `.energy`\n"
1442 " )\n"
1443 " SELECT\n"
1444 " `.energy`, lo, hi,\n" // Scale for Theta-Weights
1445 " `Signal`, `Background`/5 AS `Background`, `Excess`, `ErrExcess`, `Significance`,\n"
1446 " SQRT(`Signal`) AS ErrSignal,\n"
1447 " SQRT(`SignalW2`) AS ErrSignalW,\n"
1448 " SQRT(`Background`)/5 AS ErrBackground,\n"
1449 " SQRT(`BackgroundW2`)/5 AS ErrBackgroundW,\n"
1450 " ExcessN, ExcessW, ErrExcessN, ErrExcessW, SignalW, BackgroundW,\n"
1451 " AnalysisData.Excess/ResultMC.ExcessW*SimulatedSpectrum.CountW * 1000/(POW(10,hi)-POW(10,lo)) /Area/ObsTime / CountSim*ObsTime AS Flux,\n"
1452 " AnalysisData.Excess/ResultMC.ExcessW*SimulatedSpectrum.CountW * 1000/(POW(10,hi)-POW(10,lo)) /Area/ObsTime / CountSim*ObsTime\n"
1453 " * SQRT(\n"
1454 " + POW(AnalysisData.ErrExcess / AnalysisData.Excess, 2)\n"
1455 " + POW(ResultMC.ErrExcessW / ResultMC.ExcessW, 2)\n"
1456 " + SimulatedSpectrum.CountW2 / POW(SimulatedSpectrum.CountW,2)\n"
1457 " ) AS ErrFlux,\n"
1458 " Bias,\n"
1459 " Resolution,\n"
1460 " ResultMC.ExcessW/SimulatedSpectrum.CountW * CountSim/ObsTime AS EfficiencyW,\n"
1461 " ResultMC.ExcessN/SimulatedSpectrum.CountN AS EfficiencyN,\n"
1462 " ( POW(ResultMC.ErrExcessW/ResultMC.ExcessW, 2) + POW(SQRT(SimulatedSpectrum.CountW2)/SimulatedSpectrum.CountW, 2) )\n"
1463 " * POW(ResultMC.ExcessW/SimulatedSpectrum.CountW, 2) * CountSim/ObsTime AS ErrEfficiencyW,\n"
1464 " ( POW(ResultMC.ErrExcessN, 2) + POW(ResultMC.ExcessN, 2)/SimulatedSpectrum.CountN)/POW(SimulatedSpectrum.CountN, 2) AS ErrEfficiencyN\n"
1465 " FROM\n"
1466 " AnalysisData\n"
1467 " INNER JOIN\n"
1468 " ResultMC USING(`.energy`)\n"
1469 " INNER JOIN\n"
1470 " SimulatedSpectrum USING(`.energy`)\n"
1471 " INNER JOIN\n"
1472 " BinningEnergyEst ON `.energy`=bin\n"
1473 " CROSS JOIN\n"
1474 " ThetaSums, Area\n"
1475 " WHERE\n"
1476 " AnalysisData.Excess>0\n"
1477 " ORDER BY\n"
1478 " `.energy`\n"
1479 ")"
1480 ;
1481
1482 // [ Sa^2/a^2 + Sb^2/b^2 ] * a^2/b^2
1483 // [ (sc^2)/c^2+(sb^2)/b^2+(sa^2)/a^2 ] * a^2*b^2/c^2
1484
1485 query8.parse();
1486 for (auto it=env.cbegin(); it!=env.cend(); it++)
1487 query8.template_defaults[it->first.c_str()] = it->second.c_str();
1488
1489 //query8.template_defaults["area"] = area;
1490 //query8.template_defaults["ontime"] = resX[0]["OnTime"].data();
1491 //query8.template_defaults["count"] = resX[0]["CountN"].data();
1492
1493 if (print_queries)
1494 PrintQuery(query8.str());
1495
1496 qlog << query8 << ";\n" << endl;
1497 if (connection.connected())
1498 {
1499 cout << query8.execute().info() << endl;
1500 ShowWarnings(connection);
1501 Dump(flog, connection, "Spectrum");
1502
1503 const auto sec8 = Time().UnixTime()-start8.UnixTime();
1504 cout << "Execution time: " << sec8 << "s\n";
1505
1506
1507 const mysqlpp::StoreQueryResult res8 = connection.query("SELECT * FROM Spectrum").store();
1508
1509 if (verbose>0)
1510 {
1511 cout << " Bin Flux Error" << endl;
1512 const auto bins = binning_eest.vec();
1513 for (auto ir=res8.cbegin(); ir!=res8.cend(); ir++)
1514 {
1515 const mysqlpp::Row &row = *ir;
1516
1517 const uint32_t bin = row[".energy"];
1518 cout << setw(5) << bins[bin] << ": " << setw(10) << row["Flux"] << " " << setw(10) << row["ErrFlux"] << '\n';
1519 }
1520 cout << endl;
1521 }
1522
1523 // --------------------------------------------------------------------------
1524
1525 // Crab Nebula: 1 TeV: 3e-7 / m^2 / s / TeV
1526 // Crab Nebula: 1 TeV: 3e-11 / cm^2 / s / TeV
1527
1528 sindent indentm(mlog);
1529
1530 mlog << "void spectrum()\n";
1531 mlog << "{\n" << indent(4);
1532 mlog << "// Energy Spectrum (e.g. Crab: 3e-11 [cm^-2 s-^1 TeV^-1] @ 1TeV)\n";
1533 mlog << "TGraphErrors g;\n";
1534 mlog << "g.SetNameTitle(\"Spectrum\", \"Energy Spectrum\");\n";
1535 for (auto ir=res8.cbegin(); ir!=res8.cend(); ir++)
1536 {
1537 // This is not efficient but easier to read and safer
1538 const mysqlpp::Row &row = *ir;
1539
1540 const double hi = row["hi"];
1541 const double lo = row["lo"];
1542 const double center = (hi+lo)/2;
1543 const double flux = row["Flux"];
1544 const double error = row["ErrFlux"];
1545
1546 mlog << "g.SetPoint(g.GetN(), " << setw(8) << center << ", " << flux << ");\n";
1547 mlog << "g.SetPointError(g.GetN()-1, 0, " << error << ");\n";
1548 }
1549 mlog << "g.GetXaxis()->SetTitle(\"lg(E/TeV)\");\n";
1550 mlog << "g.GetYaxis()->SetTitle(\"dN/dE [cm^{-2} s^{-1} TeV^{-1}]\");";
1551 mlog << endl;
1552 }
1553
1554 WriteHistogram(connection, {
1555 .name = "Signal",
1556 .title = "Signal",
1557 .dir = "Eest/Measurement",
1558 .binningx = binning_eest,
1559 .binningy = {},
1560 .table = "Spectrum",
1561 .x = ".energy",
1562 .y = "",
1563 .v = "Signal",
1564 .err = "ErrSignal",
1565 .axisx = "Energy lg(E_{est}/GeV)",
1566 .axisy = "Counts",
1567 .axisz = "",
1568 .stats = true
1569 });
1570
1571 WriteHistogram(connection, {
1572 .name = "Background",
1573 .title = "Background",
1574 .dir = "Eest/Measurement",
1575 .binningx = binning_eest,
1576 .binningy = {},
1577 .table = "Spectrum",
1578 .x = ".energy",
1579 .y = "",
1580 .v = "Background",
1581 .err = "ErrBackground",
1582 .axisx = "Energy lg(E_{est}/GeV)",
1583 .axisy = "Count",
1584 .axisz = "",
1585 .stats = true
1586 });
1587
1588 WriteHistogram(connection, {
1589 .name = "Excess",
1590 .title = "Excess",
1591 .dir = "Eest/Measurement",
1592 .binningx = binning_eest,
1593 .binningy = {},
1594 .table = "Spectrum",
1595 .x = ".energy",
1596 .y = "",
1597 .v = "Excess",
1598 .err = "ErrExcess",
1599 .axisx = "Energy lg(E_{est}/GeV)",
1600 .axisy = "Signal - Background (Counts)",
1601 .axisz = "",
1602 .stats = true
1603 });
1604
1605 WriteHistogram(connection, {
1606 .name = "SignalW",
1607 .title = "SignalW",
1608 .dir = "Eest/Simulation/Weighted",
1609 .binningx = binning_eest,
1610 .binningy = {},
1611 .table = "Spectrum",
1612 .x = ".energy",
1613 .y = "",
1614 .v = "SignalW",
1615 .err = "ErrSignalW",
1616 .axisx = "Energy lg(E_{est}/GeV)",
1617 .axisy = "Weighted",
1618 .axisz = "",
1619 .stats = true
1620 });
1621
1622 WriteHistogram(connection, {
1623 .name = "BackgroundW",
1624 .title = "BackgroundW",
1625 .dir = "Eest/Simulation/Weighted",
1626 .binningx = binning_eest,
1627 .binningy = {},
1628 .table = "Spectrum",
1629 .x = ".energy",
1630 .y = "",
1631 .v = "BackgroundW",
1632 .err = "ErrBackgroundW",
1633 .axisx = "Energy lg(E_{est}/GeV)",
1634 .axisy = "Weighted",
1635 .axisz = "",
1636 .stats = true
1637 });
1638
1639 WriteHistogram(connection, {
1640 .name = "ExcessW",
1641 .title = "ExcessW",
1642 .dir = "Eest/Simulation/Weighted",
1643 .binningx = binning_eest,
1644 .binningy = {},
1645 .table = "Spectrum",
1646 .x = ".energy",
1647 .y = "",
1648 .v = "ExcessW",
1649 .err = "ErrExcessW",
1650 .axisx = "Energy lg(E_{est}/GeV)",
1651 .axisy = "Signal - Background (Weighted)",
1652 .axisz = "",
1653 .stats = true
1654 });
1655
1656 WriteHistogram(connection, {
1657 .name = "Significance",
1658 .title = "Significance",
1659 .dir = "Eest/Measurement",
1660 .binningx = binning_eest,
1661 .binningy = {},
1662 .table = "Spectrum",
1663 .x = ".energy",
1664 .y = "",
1665 .v = "Significance",
1666 .err = "",
1667 .axisx = "Energy lg(E_{est}/GeV)",
1668 .axisy = "Li/Ma Significance",
1669 .axisz = "",
1670 .stats = true
1671 });
1672
1673 WriteHistogram(connection, {
1674 .name = "Bias",
1675 .title = "Energy Bias",
1676 .dir = "Eest",
1677 .binningx = binning_eest,
1678 .binningy = {},
1679 .table = "Spectrum",
1680 .x = ".energy",
1681 .y = "",
1682 .v = "Bias",
1683 .err = "Resolution",
1684 .axisx = "Energy lg(E_{sim}/GeV)",
1685 .axisy = "lg(E_{est}/E_{sim})",
1686 .axisz = "",
1687 .stats = false
1688 });
1689
1690 WriteHistogram(connection, {
1691 .name = "Resolution",
1692 .title = "Energy Resolution",
1693 .dir = "Eest",
1694 .binningx = binning_eest,
1695 .binningy = {},
1696 .table = "Spectrum",
1697 .x = ".energy",
1698 .y = "",
1699 .v = "Resolution",
1700 .err = "",
1701 .axisx = "Energy lg(E_{sim}/GeV)",
1702 .axisy = "#sigma(lg(E_{est}/E_{sim}))",
1703 .axisz = "",
1704 .stats = false
1705 });
1706
1707 WriteHistogram(connection, {
1708 .name = "EfficiencyN",
1709 .title = "Efficiency (Counts)",
1710 .dir = "Eest/Efficiency",
1711 .binningx = binning_eest,
1712 .binningy = {},
1713 .table = "Spectrum",
1714 .x = ".energy",
1715 .y = "",
1716 .v = "EfficiencyN",
1717 .err = "ErrEfficiencyN",
1718 .axisx = "Energy lg(E_{sim}/GeV)",
1719 .axisy = "Efficiency",
1720 .axisz = "",
1721 .stats = true
1722 });
1723
1724 WriteHistogram(connection, {
1725 .name = "EfficiencyW",
1726 .title = "Efficiency (Weighted)",
1727 .dir = "Eest/Efficiency",
1728 .binningx = binning_eest,
1729 .binningy = {},
1730 .table = "Spectrum",
1731 .x = ".energy",
1732 .y = ".",
1733 .v = "EfficiencyW",
1734 .err = "ErrEfficiencyW",
1735 .axisx = "Energy lg(E_{sim}/GeV)",
1736 .axisy = "Efficiency",
1737 .axisz = "",
1738 .stats = true
1739 });
1740
1741 WriteHistogram(connection, {
1742 .name = "Spectrum",
1743 .title = "Differential Energy Spectrum",
1744 .dir = "",
1745 .binningx = binning_eest,
1746 .binningy = {},
1747 .table = "Spectrum",
1748 .x = ".energy",
1749 .y = "",
1750 .v = "Flux",
1751 .err = "ErrFlux",
1752 .axisx = "Energy lg(E/GeV)",
1753 .axisy = "dN/dE [cm^{-2} s^{-1} TeV^{-1}]",
1754 .axisz = "",
1755 .stats = false
1756 });
1757
1758
1759 // -------------------------------------------------------------------
1760 // -------------------------- 9: Threshold ---------------------------
1761
1762 cout << separator("Threshold") << '\n';
1763
1764 Time start9;
1765
1766 /* 02:calculate-threshold.sql */
1767
1768 // This query gets the analysis results versus Simulated Energy from the combined table
1769 mysqlpp::Query query9(&connection);
1770 query9 <<
1771 // SELECT SUM(CountN) AS CountN, SUM(OnTime) AS OnTime FROM ThetaHist
1772 "CREATE TEMPORARY TABLE Threshold ENGINE=Memory AS\n"
1773 "(\n"
1774 " WITH ThetaSums AS\n"
1775 " (\n"
1776 " SELECT\n"
1777 " SUM(CountN) AS CountSim,\n"
1778 " SUM(OnTime) AS ObsTime\n"
1779 " FROM\n"
1780 " ThetaHist\n"
1781 " ),\n"
1782 " Area AS\n"
1783 " (\n"
1784 " SELECT\n"
1785 " POW(MinImpactHi,2)*PI() AS Area\n"
1786 " FROM\n"
1787 " MonteCarloArea\n"
1788 " ),\n"
1789 " ResultMC AS\n" // Group AnalysisMC by EnergySim Binning
1790 " (\n"
1791 " SELECT\n"
1792 " `.energysim` AS `.energy`,\n"
1793 " ANY_VALUE(ThresholdW) AS ThresholdW,\n"
1794 " ANY_VALUE(ThresholdW2) AS ThresholdW2,\n"
1795 " ANY_VALUE(ThresholdN) AS ThresholdN,\n"
1796 " ANY_VALUE(BiasSim) AS Bias,\n"
1797 " ANY_VALUE(ResolutionSim) AS Resolution\n"
1798 " FROM\n"
1799 " AnalysisMC\n"
1800 " GROUP BY\n"
1801 " `.energy`\n"
1802 " )\n"
1803 " SELECT\n"
1804 " `.energy`, lo, hi,\n"
1805 " ThresholdW,\n"
1806 " SQRT(ThresholdW2) AS ErrThresholdW,\n"
1807 " ThresholdN,\n"
1808 " SQRT(ThresholdN) AS ErrThresholdN,\n" // Scale for Theta-Weights
1809 " ThresholdW * 1000/(POW(10,hi)-POW(10,lo)) / Area / CountSim*ObsTime AS Flux,\n"
1810 " SQRT(ThresholdW2) * 1000/(POW(10,hi)-POW(10,lo)) / Area / CountSim*ObsTime AS ErrFlux,\n"
1811 " Bias,\n"
1812 " Resolution\n"
1813 " FROM\n"
1814 " ResultMC\n"
1815 " INNER JOIN\n"
1816 " BinningEnergySim ON `.energy`=bin\n"
1817 " CROSS JOIN\n"
1818 " ThetaSums, Area\n"
1819 " WHERE\n"
1820 " ThresholdW>0 AND ThresholdW2>0\n"
1821 " ORDER BY\n"
1822 " `.energy`\n"
1823 ")";
1824
1825 query9.parse();
1826 for (auto it=env.cbegin(); it!=env.cend(); it++)
1827 query9.template_defaults[it->first.c_str()] = it->second.c_str();
1828
1829 //query9.template_defaults["area"] = area;
1830 //query9.template_defaults["ontime"] = resX[0]["OnTime"].data();
1831 //query9.template_defaults["count"] = resX[0]["CountN"].data();
1832
1833 if (print_queries)
1834 PrintQuery(query9.str());
1835
1836 qlog << query9 << ";\n" << endl;
1837 if (connection.connected())
1838 {
1839 cout << query9.execute().info() << endl;
1840 ShowWarnings(connection);
1841 Dump(flog, connection, "Threshold");
1842
1843 const auto sec9 = Time().UnixTime()-start9.UnixTime();
1844 cout << "Execution time: " << sec9 << "s\n";
1845
1846 // --------------------------------------------------------------------------
1847
1848 const mysqlpp::StoreQueryResult res9 = connection.query("SELECT * FROM Threshold").store();
1849
1850 sindent indentm(mlog, 4);
1851
1852 mlog << '\n';
1853 mlog << "// Energy Threshold\n";
1854 mlog << "TGraphErrors g;\n";
1855 mlog << "g.SetNameTitle(\"Threshold\", \"Energy Threshold\");\n";
1856 for (auto ir=res9.cbegin(); ir!=res9.cend(); ir++)
1857 {
1858 // This is not efficient but easier to read and safer
1859 const mysqlpp::Row &row = *ir;
1860
1861 const double hi = row["hi"];
1862 const double lo = row["lo"];
1863 const double center = (hi+lo)/2;
1864 const double width = pow(10, hi)-pow(10, lo);
1865 const double flux = row["Flux"] /width;
1866 const double error = row["ErrFlux"]/width;
1867
1868 mlog << "g.SetPoint(g.GetN(), " << setw(8) << center << ", " << flux << ");\n";
1869 mlog << "g.SetPointError(g.GetN()-1, 0, " << error << ");\n";
1870 }
1871 mlog << "g.GetXaxis()->SetTitle(\"lg(E/TeV)\");\n";
1872 mlog << "g.GetYaxis()->SetTitle(\"dN/dE [cm^{-2} s^{-1} TeV^{-1}]\");\n";
1873 mlog << indent(0) << "}" << endl;
1874 }
1875
1876 WriteHistogram(connection, {
1877 .name = "SimExcessW",
1878 .title = "Weighted Simulated Excess",
1879 .dir = "Esim/Simulation/Weighted",
1880 .binningx = binning_esim,
1881 .binningy = {},
1882 .table = "Threshold",
1883 .x = ".energy",
1884 .y = "",
1885 .v = "ThresholdW",
1886 .err = "ErrThresholdW",
1887 .axisx = "Energy lg(E_{sim}/GeV)",
1888 .axisy = "Weighted Counts",
1889 .axisz = "",
1890 .stats = true
1891 });
1892
1893 WriteHistogram(connection, {
1894 .name = "SimExcessN",
1895 .title = "Simulated Excess",
1896 .dir = "Esim/Simulation/Counts",
1897 .binningx = binning_esim,
1898 .binningy = {},
1899 .table = "Threshold",
1900 .x = ".energy",
1901 .y = "",
1902 .v = "ThresholdN",
1903 .err = "ErrThresholdN",
1904 .axisx = "Energy lg(E_{sim}/GeV)",
1905 .axisy = "Counts",
1906 .axisz = "",
1907 .stats = true
1908 });
1909
1910 WriteHistogram(connection, {
1911 .name = "SimSpectrumW",
1912 .title = "Weighted Simulated Excess Spectrum",
1913 .dir = "Esim/Simulation/Weighted",
1914 .binningx = binning_esim,
1915 .binningy = {},
1916 .table = "Threshold",
1917 .x = ".energy",
1918 .y = "",
1919 .v = "Flux",
1920 .err = "ErrFlux",
1921 .axisx = "Energy lg(E_{sim}/GeV)",
1922 .axisy = "dN/dE [cm^{-2} s^{-1} TeV^{-1}]",
1923 .axisz = "",
1924 .stats = true
1925 });
1926
1927 WriteHistogram(connection, {
1928 .name = "Bias",
1929 .title = "Energy Bias",
1930 .dir = "Esim",
1931 .binningx = binning_esim,
1932 .binningy = {},
1933 .table = "Threshold",
1934 .x = ".energy",
1935 .y = "",
1936 .v = "Bias",
1937 .err = "Resolution",
1938 .axisx = "Energy lg(E_{sim}/GeV)",
1939 .axisy = "lg(E_{est}/E_{sim})",
1940 .axisz = "",
1941 .stats = false
1942 });
1943
1944 WriteHistogram(connection, {
1945 .name = "Resolution",
1946 .title = "Energy Resolution",
1947 .dir = "Esim",
1948 .binningx = binning_esim,
1949 .binningy = {},
1950 .table = "Threshold",
1951 .x = ".energy",
1952 .y = "",
1953 .v = "Resolution",
1954 .err = "",
1955 .axisx = "Energy lg(E_{sim}/GeV)",
1956 .axisy = "#sigma(lg(E_{est}/E_{sim}))",
1957 .axisz = "",
1958 .stats = false
1959 });
1960
1961
1962 // -------------------------------------------------------------------
1963 // -------------------------- 10: Migration --------------------------
1964
1965 cout << separator("Migration") << '\n';
1966
1967 Time start10;
1968
1969 /* 02:obtain-migration-matrix.sql */
1970
1971 // This query gets the analysis results versus Estimated Energy from the combined table
1972 mysqlpp::Query query10(&connection);
1973 query10 <<
1974 "CREATE TEMPORARY TABLE Migration ENGINE=Memory AS\n"
1975 "(\n"
1976 " SELECT\n"
1977 " `.energyest`,\n"
1978 " `.energysim`,\n"
1979 " BinningEnergySim.lo AS EsimLo,\n"
1980 " BinningEnergySim.hi AS EsimHi,\n"
1981 " BinningEnergyEst.lo AS EestLo,\n"
1982 " BinningEnergyEst.hi AS EestHi,\n"
1983 " ANY_VALUE(MigrationW) AS MigrationW,\n"
1984 " ANY_VALUE(MigrationN) AS MigrationN\n"
1985 // FIXME: Errors
1986 " FROM\n"
1987 " AnalysisMC\n"
1988 " INNER JOIN\n"
1989 " BinningEnergyEst ON `.energyest`=BinningEnergyEst.bin\n"
1990 " INNER JOIN\n"
1991 " BinningEnergySim ON `.energysim`=BinningEnergySim.bin\n"
1992 " GROUP BY\n"
1993 " `.energyest`, `.energysim`\n"
1994 " ORDER BY\n"
1995 " `.energyest`, `.energysim`\n"
1996 ")";
1997
1998 if (print_queries)
1999 PrintQuery(query10.str());
2000
2001 qlog << query10 << ";\n" << endl;
2002 if (connection.connected())
2003 {
2004 cout << query10.execute().info() << endl;
2005 ShowWarnings(connection);
2006 Dump(flog, connection, "Migration");
2007
2008 const auto sec10 = Time().UnixTime()-start10.UnixTime();
2009 cout << "Execution time: " << sec10 << "s\n\n";
2010 }
2011
2012 WriteHistogram(connection, {
2013 .name = "MigrationN",
2014 .title = "Energy Migration",
2015 .dir = "",
2016 .binningx = binning_esim,
2017 .binningy = binning_eest,
2018 .table = "Migration",
2019 .x = ".energysim",
2020 .y = ".energyest",
2021 .v = "MigrationN",
2022 .err = "",
2023 .axisx = "Energy lg(E_{sim}/GeV)",
2024 .axisy = "Energy lg(E_{est}/GeV)",
2025 .axisz = "Counts",
2026 .stats = false
2027 });
2028
2029 WriteHistogram(connection, {
2030 .name = "MigrationW",
2031 .title = "Energy Migration",
2032 .dir = "",
2033 .binningx = binning_esim,
2034 .binningy = binning_eest,
2035 .table = "Migration",
2036 .x = ".energysim",
2037 .y = ".energyest",
2038 .v = "MigrationW",
2039 .err = "",
2040 .axisx = "Energy lg(E_{sim}/GeV)",
2041 .axisy = "Energy lg(E_{est}/GeV)",
2042 .axisz = "Weighted Counts",
2043 .stats = false
2044 });
2045
2046
2047
2048 // -------------------------------------------------------------------
2049 // --------------------------- 11: Summary ---------------------------
2050
2051 cout << separator("Summary") << '\n';
2052 const auto sec = Time().UnixTime()-start.UnixTime();
2053 cout << "Total execution time: " << sec << "s\n" << endl;
2054
2055 return 0;
2056}
Note: See TracBrowser for help on using the repository browser.