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

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