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

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