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

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