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

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