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

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