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

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