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

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