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

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