Changeset 19895
- Timestamp:
- 12/14/19 19:39:38 (5 years ago)
- File:
-
- 1 edited
-
trunk/FACT++/src/spectrum.cc (modified) (32 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/FACT++/src/spectrum.cc
r19889 r19895 17 17 #include "TH1.h" 18 18 #include "TH2.h" 19 #include "TRolke.h" 20 #include "TFeldmanCousins.h" 19 21 #endif 20 22 … … 130 132 #endif 131 133 , "Database link as in\n\tuser:password@server[:port]/database[?compress=0|1].") 132 ("out,o", var<string>(conf.GetName()), "Defines the prefix (with path) of the output files.") 133 ("dry-run", po_switch(), "Only write the queries to the query.log file. Internally overwrites uri with an empty string.") 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).") 134 137 ; 135 138 136 139 po::options_description binnings("Binnings"); 137 140 binnings.add_options() 138 ("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)")139 ("theta-bin", vars<double>(), "Add a bin-edge to the theta binning (degree)")140 ("e sim", var<Binning>(Binning(15, 2, 5)), "Add equidistant bins in log10 simulated energy. Syntax: N,lo,hi (Number N of equidistant bins between lo and hi)")141 ("e sim-bin", vars<double>(), "Add a bin-edge to the binnnig in log10 simulated enegry")142 ("e est",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)")143 ("e est-bin",vars<double>(), "Add a bin-edge to the binning in log10 estimated enegry")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") 144 147 ; 145 148 … … 147 150 queries.add_options() 148 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.") 149 ("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.")150 ("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.")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.") 151 154 ; 152 155 … … 162 165 po::options_description debug("Debug options"); 163 166 debug.add_options() 164 ("print-connection", po_switch(), "Print database connection information") 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") 165 169 #ifdef HAVE_HIGHLIGHT 166 ("print-queries", po_ switch(),"Print all queries to the console. They are automatically highlighted.")170 ("print-queries", po_bool(), "Print all queries to the console. They are automatically highlighted.") 167 171 #else 168 ("print-queries", po_ switch(),"Print all queries to the console. (For highlighting recompile with 'libsource-highlight-dev' installed)")172 ("print-queries", po_bool(), "Print all queries to the console. (For highlighting recompile with 'libsource-highlight-dev' installed)") 169 173 #endif 174 ("mc-only", po_bool(), "Do not run an data realated queries (except observation times)") 170 175 ("verbose,v", var<uint16_t>(1), "Verbosity (0: quiet, 1: default, 2: more, 3, ...)") 171 176 ; … … 327 332 srchilite::SourceHighlight sourceHighlight("esc256.outlang"); 328 333 sourceHighlight.setStyleFile("esc256.style"); 334 sourceHighlight.setGenerateLineNumbers(); 335 sourceHighlight.setLineNumberDigits(3); 336 //sourceHighlight.setLineNumberPad(' ') 329 337 sourceHighlight.highlight(qstr, cout, "sql.lang"); 330 338 cout << endl; … … 341 349 "(\n" 342 350 " bin SMALLINT UNSIGNED NOT NULL,\n" 343 " lo FLOATNOT NULL,\n"344 " hi FLOATNOT NULL,\n"351 " lo DOUBLE NOT NULL,\n" 352 " hi DOUBLE NOT NULL,\n" 345 353 " PRIMARY KEY (bin) USING HASH\n" 346 354 ")"; … … 349 357 if (connection.connected()) 350 358 query0.execute(); 351 352 //connection.query("ALTER TABLE Binning"+name+" AUTO_INCREMENT=0").execute();353 359 354 360 mysqlpp::Query query1(&connection); … … 426 432 }; 427 433 434 #ifdef HAVE_ROOT 435 436 TH1 *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 480 void 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 428 495 void WriteHistogram(Database &connection, const Histogram &hist) 429 496 { … … 431 498 if (!connection.connected()) 432 499 return; 500 501 TH1 *h = CreateHistogram(hist); 433 502 434 503 mysqlpp::Query query(&connection); … … 449 518 const mysqlpp::StoreQueryResult res = query.store(); 450 519 451 const auto vecx = hist.binningx.vec();452 const auto vecy = hist.binningy.vec();453 454 TH1 *h = 0;455 456 if (hist.y.empty())457 {458 h = hist.binningx.equidist ?459 new TH1D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.front(), vecx.back()) :460 new TH1D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.data());461 }462 else463 {464 if (hist.binningy.equidist)465 {466 h = hist.binningx.equidist ?467 new TH2D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.front(), vecx.back(), vecy.size()-1, vecy.front(), vecy.back()) :468 new TH2D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.data(), vecy.size()-1, vecy.front(), vecy.back());469 }470 else471 {472 h = hist.binningx.equidist ?473 new TH2D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.front(), vecx.back(), vecy.size()-1, vecy.front(), vecy.back()) :474 new TH2D(hist.name.c_str(), hist.title.c_str(), vecx.size()-1, vecx.data(), vecy.size()-1, vecy.front(), vecy.back());475 }476 }477 478 h->SetXTitle(hist.axisx.c_str());479 h->SetYTitle(hist.axisy.c_str());480 h->SetZTitle(hist.axisz.c_str());481 482 h->SetMarkerStyle(kFullDotMedium);483 484 h->SetStats(hist.stats);485 486 520 for (auto ir=res.cbegin(); ir!=res.cend(); ir++) 487 521 { 488 522 const auto &row = *ir; 489 490 const uint32_t x = row["X"];491 const double v = row["V"];492 if (x>uint32_t(h->GetNbinsX()+1))493 continue;494 523 495 524 try 496 525 { 497 const uint32_t y = row["Y"]; 498 if (y>uint32_t(h->GetNbinsY()+1)) 526 const uint32_t x = row["X"]; 527 const double v = row["V"]; 528 if (x>uint32_t(h->GetNbinsX()+1)) 499 529 continue; 500 530 501 h->SetBinContent(x, y, v);502 503 }504 catch (const mysqlpp::BadFieldName &)505 {506 h->SetBinContent(x, v);507 531 try 508 532 { 509 h->SetBinError(x, double(row["E"])); 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 510 539 } 511 540 catch (const mysqlpp::BadFieldName &) 512 541 { 542 h->SetBinContent(x, v); 543 try 544 { 545 h->SetBinError(x, double(row["E"])); 546 } 547 catch (const mysqlpp::BadFieldName &) 548 { 549 } 513 550 } 514 551 } 515 } 516 517 gFile->cd(); 518 TDirectory *dir = gFile->GetDirectory(hist.dir.c_str()); 519 if (dir) 520 dir->cd(); 521 else 522 { 523 gFile->mkdir(hist.dir.c_str()); 524 gFile->cd(hist.dir.c_str()); 525 } 526 h->Write(); 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 563 void 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); 527 572 delete h; 528 573 #endif … … 581 626 // ----------------------------- Evaluate options -------------------------- 582 627 583 const string uri = conf.Get<bool>("dry-run") ? "" : conf.Get<string>("uri"); 584 const string out = conf.Get<string>("out"); 585 const uint16_t verbose = conf.Get<uint16_t>("verbose"); 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"); 586 633 587 634 const bool print_connection = conf.Get<bool>("print-connection"); 588 635 const bool print_queries = conf.Get<bool>("print-queries"); 589 590 Binning binning_theta = conf.Get<Binning>("theta")+conf.Vec<double>("theta-bin"); 591 Binning binning_esim = conf.Get<Binning>("esim"); 592 Binning binning_eest = conf.Get<Binning>("eest"); 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"); 593 641 594 642 cout << '\n'; 595 cout << "Binning 'theta': " << binning_theta.str() << endl;596 cout << "Binning ' Esim': " << binning_esim.str() << endl;597 cout << "Binning ' Eest': " << binning_eest.str() << endl;643 cout << "Binning 'theta': " << binning_theta.str() << endl; 644 cout << "Binning 'dense': " << binning_dense.str() << endl; 645 cout << "Binning 'sparse': " << binning_sparse.str() << endl; 598 646 599 647 const uint16_t source_key = conf.Get<uint16_t>("source-key"); … … 604 652 605 653 cout << "\n"; 606 const string analysis_sql = CreateResource(conf, "analysis", "analysis.sql", RESOURCE(std::string, spectrum_analysis_sql)); 607 const string data_sql = CreateResource(conf, "data", "data.sql", RESOURCE(std::string, spectrum_data_sql)); 608 const string simulation_sql = CreateResource(conf, "simulation", "simulation.sql", RESOURCE(std::string, spectrum_simulation_sql)); 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); 609 660 cout << endl; 610 661 611 const string str zd= binning_theta.list();612 const string str e = binning_eest.list();613 const string str th = binning_esim.list();662 const string str_theta = binning_theta.list(); 663 const string str_dense = binning_dense.list(); 664 const string str_sparse = binning_sparse.list(); 614 665 615 666 // ------------------------------------------------------------------------- … … 696 747 cout << separator("Binnings") << '\n'; 697 748 698 CreateBinning(connection, qlog, binning_theta, "Theta");699 CreateBinning(connection, qlog, binning_ eest, "EnergyEst");700 CreateBinning(connection, qlog, binning_ esim, "EnergySim");749 CreateBinning(connection, qlog, binning_theta, "Theta"); 750 CreateBinning(connection, qlog, binning_dense, "Energy_dense"); 751 CreateBinning(connection, qlog, binning_sparse, "Energy_sparse"); 701 752 702 753 Dump(flog, connection, "BinningTheta"); 703 Dump(flog, connection, "BinningEnergy Est");704 Dump(flog, connection, "BinningEnergy Sim");754 Dump(flog, connection, "BinningEnergy_dense"); 755 Dump(flog, connection, "BinningEnergy_sparse"); 705 756 706 757 // ------------------------------------------------------------------- … … 719 770 " FileId INT UNSIGNED NOT NULL,\n" 720 771 " PRIMARY KEY (FileId) USING HASH\n" 721 ") ENGINE=M emory\n"772 ") ENGINE=MEMORY\n" 722 773 "AS\n" 723 774 "(\n" … … 770 821 "(\n" 771 822 " `.theta` SMALLINT UNSIGNED NOT NULL,\n" 772 " OnTime FLOATNOT NULL,\n"823 " OnTime DOUBLE NOT NULL,\n" 773 824 " PRIMARY KEY (`.theta`) USING HASH\n" 774 ") ENGINE=M emory\n"825 ") ENGINE=MEMORY\n" 775 826 "AS\n" 776 827 "(\n" … … 792 843 query2.template_defaults[it->first.c_str()] = it->second.c_str(); 793 844 794 query2.template_defaults["bins"] = str zd.c_str();845 query2.template_defaults["bins"] = str_theta.c_str(); 795 846 796 847 if (print_queries) … … 808 859 } 809 860 810 /*811 812 SELECT813 min(first_number) as first_number,814 max(last_number) as last_number815 from816 (817 select818 first_number,819 last_number820 sum(grp) over(order by first_number) as grp821 from822 (823 select824 bin AS first_number,825 LEAD(bin) AS last_number826 IF(first_number <> lag(bin, 2) over (order by bin) + 1, 1, 0) as grp827 from828 t1829 )830 )831 group by grp832 order by 1833 */834 835 /*836 WITH Hallo AS (WITH MyTest AS (WITH MyTable AS (SELECT bin FROM Test837 WHERE cnt>0 GROUP BY bin ORDER BY bin) SELECT bin AS first_number,838 LAG(bin) OVER(ORDER BY bin) AS prev, LEAD(bin) OVER(ORDER BY bin) AS839 last_number, IF(bin <> (LAG(bin) OVER(ORDER BY bin)+1), 1, 0) AS840 grp_prev, IF(bin <> (LEAD(bin) OVER(ORDER BY bin)-1), 1, 0) as grp_nxt841 FROM MyTable ) SELECT first_number, last_number, sum(grp_nxt)842 over(order by first_number) as grp FROM MyTest) SELECT * FROM Hallo843 844 845 WITH MyTable AS846 (SELECT bin FROM Test WHERE cnt>0 GROUP BY bin ORDER BY bin)847 SELECT bin, LAG(bin) OVER(ORDER BY bin) AS prev, LEAD(bin) OVER(ORDER BY bin) AS nxt,848 IF(bin = (LAG(bin) OVER(ORDER BY bin)+1), 1, 0) AS grp_prev, IF(bin = (LEAD(bin) OVER(ORDER BY bin)-1), 1, 0) as grp_nxt FROM MyTable849 850 "-1","0"851 "0","0"852 "2","3" 2-4853 "3","5"854 "4","6"855 "5","0"856 "6","0"857 "7","7" 7-7858 "8","0"859 "9","9" 9-12860 "10","10"861 "11","11"862 "12","12"863 "13","0"864 "14","0"865 "15","0"866 "16","16" 16-19867 "17","17"868 "18","18"869 "19","19"870 */871 872 861 // ------------------------------------------------------------------- 873 862 // -------------------------- MonteCarloFiles ------------------------ … … 878 867 Time start3; 879 868 880 /* 02:get-monte-carlo-files.sql */881 869 mysqlpp::Query query3(&connection); 882 870 query3 << … … 885 873 " FileId INT UNSIGNED NOT NULL,\n" 886 874 " PRIMARY KEY (FileId) USING HASH\n" 887 ") ENGINE=M emory\n"875 ") ENGINE=MEMORY\n" 888 876 "AS\n" 889 877 "(\n" … … 897 885 " factmc.RunInfoMC\n" 898 886 " ON\n" 899 //" ThetaMin BETWEEN lo AND hi OR ThetaMax BETWEEN lo AND hi\n" // Includes BOTH edges900 887 " (ThetaMin>=lo AND ThetaMin<hi) OR (ThetaMax>lo AND ThetaMax<=hi)\n" 901 888 " WHERE\n" … … 905 892 " FileId\n" // In order: faster 906 893 ")"; 907 /*908 " SELECT\n"909 " FileId\n"910 " FROM\n"911 " factmc.RunInfoMC\n"912 " WHERE\n"913 " PartId=1 AND\n"914 " FileId%%2=0 AND\n"915 " (%0:range)\n"916 ")";917 */918 894 919 895 query3.parse(); … … 941 917 cout << separator("MonteCarloArea") << '\n'; 942 918 943 Time start3b; 944 945 /* 02:get-monte-carlo-files.sql */ 946 mysqlpp::Query query3b(&connection); 947 query3b << 948 "CREATE TEMPORARY TABLE MonteCarloArea ENGINE=Memory\n" 919 Time start4; 920 921 mysqlpp::Query query4(&connection); 922 query4 << 923 "CREATE TEMPORARY TABLE MonteCarloArea ENGINE=MEMORY\n" 949 924 "AS\n" 950 925 "(\n" … … 964 939 ")"; 965 940 966 query 3b.parse();941 query4.parse(); 967 942 for (auto it=env.cbegin(); it!=env.cend(); it++) 968 query 3b.template_defaults[it->first.c_str()] = it->second.c_str();943 query4.template_defaults[it->first.c_str()] = it->second.c_str(); 969 944 970 945 if (print_queries) 971 PrintQuery(query 3b.str());972 973 qlog << query 3b<< ";\n" << endl;946 PrintQuery(query4.str()); 947 948 qlog << query4 << ";\n" << endl; 974 949 if (connection.connected()) 975 950 { 976 cout << query 3b.execute().info() << endl;951 cout << query4.execute().info() << endl; 977 952 ShowWarnings(connection); 978 953 if (Dump(flog, connection, "MonteCarloArea")!=1) … … 982 957 } 983 958 984 const auto sec 3b = Time().UnixTime()-start3b.UnixTime();985 cout << "Execution time: " << sec 3b<< "s\n\n";986 } 987 988 989 // ------------------------------------------------------------------- 990 // ------------------------ ThetaDistribution------------------------991 // ------------------------------------------------------------------- 992 993 cout << separator(" ThetaDistribution") << '\n';994 995 Time start 4;996 997 / * 02:get-theta-distribution.sql */998 mysqlpp::Query query 4(&connection);999 query 4<<1000 "CREATE TEMPORARY TABLE ThetaDistribution\n"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" 1001 976 "(\n" 1002 " `.theta` SMALLINT UNSIGNED NOT NULL,\n"1003 " lo DOUBLE NOT NULL COMMENT 'Lower edge of zenith distance bin in degree',\n"1004 " hi DOUBLE NOT NULL COMMENT 'Upper edge of zenith distance bin in degree',\n"1005 " CountN INT UNSIGNEDNOT NULL,\n"1006 " ErrCountNDOUBLE NOT NULL,\n"1007 " OnTime FLOATNOT NULL,\n"1008 " ZdWeight DOUBLE NOT NULL COMMENT 'tau(delta theta)',\n"1009 " ErrZdWeight DOUBLE NOT NULL COMMENT 'sigma(tau)',\n"1010 " PRIMARY KEY (`.theta`)USING HASH\n"1011 ") ENGINE=M emory\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" 1012 987 "AS\n" 1013 988 "(\n" 1014 " WITH EventCountAS\n"989 " WITH Table0 AS\n" 1015 990 " (\n" 1016 991 " SELECT\n" 1017 " INTERVAL(DEGREES(Theta), %100:bins) AS `.theta`,\n" 1018 " COUNT(*) AS CountN\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" 1019 996 " FROM\n" 1020 997 " MonteCarloFiles\n" 1021 998 " LEFT JOIN\n" 1022 " factmc.OriginalMC USING(FileId)\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" 1023 1084 " GROUP BY\n" 1024 1085 " `.theta`\n" 1025 1086 " )\n" 1026 1087 " SELECT\n" 1027 " `.theta`, lo, hi,\n"1088 " `.theta`,\n" 1028 1089 " CountN,\n" 1029 1090 " SQRT(CountN) AS ErrCountN,\n" … … 1034 1095 " ObservationTime\n" 1035 1096 " LEFT JOIN\n" 1036 " EventCount USING(`.theta`)\n"1097 " ThetaCount USING (`.theta`)\n" 1037 1098 " LEFT JOIN\n" 1038 1099 " BinningTheta ON `.theta`=bin\n" … … 1041 1102 ")"; 1042 1103 1043 query 4.parse();1104 query6.parse(); 1044 1105 for (auto it=env.cbegin(); it!=env.cend(); it++) 1045 query4.template_defaults[it->first.c_str()] = it->second.c_str(); 1046 1047 query4.template_defaults["bins"] = strzd.c_str(); 1106 query6.template_defaults[it->first.c_str()] = it->second.c_str(); 1048 1107 1049 1108 if (print_queries) 1050 PrintQuery(query 4.str());1051 1052 qlog << query 4<< ";\n" << endl;1109 PrintQuery(query6.str()); 1110 1111 qlog << query6 << ";\n" << endl; 1053 1112 if (connection.connected()) 1054 1113 { 1055 cout << query 4.execute().info() << endl;1114 cout << query6.execute().info() << endl; 1056 1115 ShowWarnings(connection); 1057 Dump(flog, connection, "ThetaDistribution"); 1058 1059 const auto sec4 = Time().UnixTime()-start4.UnixTime(); 1060 cout << "Execution time: " << sec4 << "s\n"; 1061 1062 1063 if (verbose>0) 1064 { 1065 const mysqlpp::StoreQueryResult res4 = connection.query("SELECT * FROM ThetaDistribution").store(); 1066 1067 cout << " Bin MC Counts OnTime ZdWeight\n"; 1068 const auto bins = binning_theta.vec(); 1069 for (auto ir=res4.cbegin(); ir!=res4.cend(); ir++) 1070 { 1071 const mysqlpp::Row &row = *ir; 1072 1073 const uint32_t bin = row[".theta"]; 1074 cout << setw(5) << bins[bin] << ": " << setw(10) << row["CountN"] << " " << setw(10) << row["OnTime"] << " " << setw(10) << row["ZdWeight"] << '\n'; 1075 } 1076 cout << endl; 1077 } 1116 Dump(flog, connection, "ThetaDist"); 1117 1118 const auto sec6 = Time().UnixTime()-start6.UnixTime(); 1119 cout << "Execution time: " << sec6 << "s\n\n"; 1078 1120 } 1079 1121 … … 1084 1126 .binningx = binning_theta, 1085 1127 .binningy = {}, 1086 .table = "ThetaDist ribution",1128 .table = "ThetaDist", 1087 1129 .x = ".theta", 1088 1130 .y = "", … … 1101 1143 .binningx = binning_theta, 1102 1144 .binningy = {}, 1103 .table = "ThetaDist ribution",1145 .table = "ThetaDist", 1104 1146 .x = ".theta", 1105 1147 .y = "", … … 1118 1160 .binningx = binning_theta, 1119 1161 .binningy = {}, 1120 .table = "ThetaDist ribution",1162 .table = "ThetaDist", 1121 1163 .x = ".theta", 1122 1164 .y = "", … … 1129 1171 }); 1130 1172 1131 1132 // ------------------------------------------------------------------- 1133 // ------------------------- 5: AnalysisData ------------------------- 1134 // ------------------------------------------------------------------- 1135 1136 cout << separator("AnalysisData") << '\n'; 1137 1138 Time start5; 1139 1140 /* 02:analyze-data.sql */ 1141 mysqlpp::Query query5(&connection); 1142 sindent indent5(query5); 1143 query5 << 1144 "CREATE TEMPORARY TABLE AnalysisData\n" 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" 1145 1195 "(\n" 1146 " `.energy` SMALLINT UNSIGNED NOT NULL,\n" 1147 " `Signal` DOUBLE NOT NULL,\n" 1148 " `Background` DOUBLE NOT NULL,\n" 1149 " `Excess` DOUBLE NOT NULL,\n" 1150 " `Significance` DOUBLE NOT NULL,\n" 1151 " `ErrExcess` DOUBLE NOT NULL,\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" 1152 1353 " PRIMARY KEY (`.energy`) USING HASH\n" 1153 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" : "1000/(POW(10,hi)-POW(10,lo))"; 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"*/ 1154 1724 "AS\n" 1155 1725 "(\n" … … 1160 1730 " Result AS\n" 1161 1731 " (\n" << indent(6) 1162 << ifstream(data_sql).rdbuf() << indent(0)<< // Must end with EOL and not in the middle of a comment1732 << data_sql << indent(0) << // Must end with EOL and not in the middle of a comment 1163 1733 " )\n" 1164 1734 " SELECT * FROM Result\n" 1165 1735 ")"; 1166 1736 1167 query 5.parse();1737 query12.parse(); 1168 1738 for (auto it=env.cbegin(); it!=env.cend(); it++) 1169 query 5.template_defaults[it->first.c_str()] = it->second.c_str();1739 query12.template_defaults[it->first.c_str()] = it->second.c_str(); 1170 1740 1171 1741 //query5.template_defaults["columns"] = "FileId, EvtNumber,"; 1172 query5.template_defaults["columns"] = ""; 1173 query5.template_defaults["files"] = "DataFiles"; 1174 query5.template_defaults["runinfo"] = "factdata.RunInfo"; 1175 query5.template_defaults["events"] = "factdata.Images"; 1176 query5.template_defaults["positions"] = "factdata.Position"; 1177 1178 query5.template_defaults["bins"] = stre.c_str(); 1179 query5.template_defaults["estimator"] = estimator.c_str(); 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(); 1180 1751 1181 1752 if (print_queries) 1182 PrintQuery(query 5.str());1183 1184 qlog << query 5<< ";\n" << endl;1753 PrintQuery(query12.str()); 1754 1755 qlog << query12 << ";\n" << endl; 1185 1756 if (connection.connected()) 1186 1757 { 1187 cout << query 5.execute().info() << endl;1758 cout << query12.execute().info() << endl; 1188 1759 ShowWarnings(connection); 1189 Dump(flog, connection, "AnalysisData"); 1190 1191 const auto sec5 = Time().UnixTime()-start5.UnixTime(); 1192 cout << "Execution time: " << sec5 << "s\n"; 1193 1194 if (verbose>0) 1195 { 1196 const mysqlpp::StoreQueryResult res5 = connection.query("SELECT * FROM AnalysisData").store(); 1197 1198 cout << " Bin Signal Background Excess Significance Error" << endl; 1199 for (auto row=res5.cbegin(); row!=res5.cend(); row++) 1200 { 1201 for (auto it=row->begin(); it!=row->end(); it++) 1202 cout << setw(10) << *it << " "; 1203 cout << '\n'; 1204 } 1205 cout << endl; 1206 } 1207 } 1208 1209 1210 // ------------------------------------------------------------------- 1211 // ----------------------------- Triggered --------------------------- 1212 // ------------------------------------------------------------------- 1213 1214 cout << separator("Triggered") << '\n'; 1215 1216 Time start6a; 1217 1218 /* 02:analyze-simulation.sql */ 1219 1220 // This table combines the analysis results vs. Binning in Estimated Energy and Simulated Energy 1221 mysqlpp::Query query6a(&connection); 1222 query6a << 1223 "CREATE TEMPORARY TABLE Triggered\n" 1224 "(\n" 1225 " `.energy` SMALLINT UNSIGNED NOT NULL,\n" 1226 " CountN DOUBLE NOT NULL,\n" 1227 " CountW DOUBLE NOT NULL,\n" 1228 " CountW2 DOUBLE NOT NULL,\n" 1229 " INDEX (`.energy`) USING HASH\n" 1230 ") ENGINE=Memory\n" 1231 "AS\n" 1232 "(\n" 1233 " WITH TriggeredEvents AS\n" 1234 " (\n" 1235 " SELECT\n" 1236 " INTERVAL(Zd, %100:theta) AS `.theta`,\n" 1237 " INTERVAL(log10(Energy), %101:energy) AS `.energy`,\n" 1238 " (%102:spectrum)/pow(Energy, SpectralIndex) AS SpectralWeight\n" 1239 " FROM\n" 1240 " MonteCarloFiles\n" 1241 " LEFT JOIN\n" 1242 " factmc.RunInfoMC USING (FileId)\n" 1243 " LEFT JOIN\n" 1244 " factmc.EventsMC USING (FileId)\n" 1245 " )\n" 1246 " SELECT\n" 1247 " `.energy`,\n" 1248 " COUNT(*) AS CountN,\n" 1249 " SUM(ZdWeight*SpectralWeight) AS CountW,\n" 1250 " SUM(POW(ErrZdWeight*SpectralWeight,2)) AS `CountW2`\n" 1251 " FROM\n" 1252 " TriggeredEvents\n" 1253 " INNER JOIN\n" 1254 " ThetaDistribution USING(`.theta`)\n" 1255 " GROUP BY\n" 1256 " `.energy`\n" 1257 ")"; 1258 1259 query6a.parse(); 1260 for (auto it=env.cbegin(); it!=env.cend(); it++) 1261 query6a.template_defaults[it->first.c_str()] = it->second.c_str(); 1262 1263 query6a.template_defaults["energy"] = strth.c_str(); 1264 query6a.template_defaults["theta"] = strzd.c_str(); 1265 query6a.template_defaults["spectrum"] = spectrum.c_str(); 1266 1267 if (print_queries) 1268 PrintQuery(query6a.str()); 1269 1270 qlog << query6a << ";\n" << endl; 1271 if (connection.connected()) 1272 { 1273 cout << query6a.execute().info() << endl; 1274 ShowWarnings(connection); 1275 Dump(flog, connection, "Triggered"); 1276 1277 const auto sec6a = Time().UnixTime()-start6a.UnixTime(); 1278 cout << "Execution time: " << sec6a << "s\n\n"; 1279 } 1280 1281 // ------------------------------------------------------------------- 1282 // ----------------------------- ResultMC ---------------------------- 1283 // ------------------------------------------------------------------- 1284 1285 cout << separator("ResultMC") << '\n'; 1286 1287 Time start6; 1288 1289 /* 02:analyze-simulation.sql */ 1290 1291 // This table combines the analysis results vs. Binning in Estimated Energy and Simulated Energy 1292 mysqlpp::Query query6(&connection); 1293 sindent indent6(query6); 1294 query6 << 1295 "CREATE TEMPORARY TABLE AnalysisMC\n" 1296 "(\n" 1297 " `.energyest` SMALLINT UNSIGNED NOT NULL,\n" 1298 " `.energysim` SMALLINT UNSIGNED NOT NULL,\n" 1299 " SignalW DOUBLE NOT NULL,\n" 1300 " BackgroundW DOUBLE NOT NULL,\n" 1301 " ThresholdW DOUBLE NOT NULL,\n" 1302 " SignalN DOUBLE NOT NULL,\n" 1303 " BackgroundN DOUBLE NOT NULL,\n" 1304 " ThresholdN DOUBLE NOT NULL,\n" 1305 " ExcessW DOUBLE NOT NULL,\n" 1306 " ExcessN DOUBLE NOT NULL,\n" 1307 " ErrExcessW DOUBLE NOT NULL,\n" 1308 " ErrExcessN DOUBLE NOT NULL,\n" 1309 " BiasEst DOUBLE NOT NULL,\n" 1310 " BiasSim DOUBLE NOT NULL,\n" 1311 " ResolutionEst DOUBLE,\n" 1312 " ResolutionSim DOUBLE,\n" 1313 " INDEX (`.energyest`) USING HASH,\n" 1314 " INDEX (`.energysim`) USING HASH\n" 1315 ") ENGINE=Memory\n" 1316 "AS\n" 1317 "(\n" 1318 " WITH Excess AS\n" 1319 " (\n" << indent(6) 1320 << ifstream(analysis_sql).rdbuf() << indent(0) << 1321 " ),\n" << 1322 " Result AS\n" 1323 " (\n" << indent(6) 1324 << ifstream(simulation_sql).rdbuf() << indent(0) << // Must end with EOL and not in the middle of a comment 1325 " )\n" 1326 " SELECT * FROM Result\n" 1327 ")"; 1328 1329 query6.parse(); 1330 for (auto it=env.cbegin(); it!=env.cend(); it++) 1331 query6.template_defaults[it->first.c_str()] = it->second.c_str(); 1332 1333 //query6.template_defaults["columns"] = "FileId, EvtNumber, CorsikaNumReuse,"; 1334 query6.template_defaults["columns"] = "Zd, Energy, SpectralIndex,"; 1335 query6.template_defaults["files"] = "MonteCarloFiles"; 1336 query6.template_defaults["runinfo"] = "factmc.RunInfoMC"; 1337 query6.template_defaults["events"] = "factmc.EventsMC"; 1338 query6.template_defaults["positions"] = "factmc.PositionMC"; 1339 1340 query6.template_defaults["energyest"] = stre.c_str(); 1341 query6.template_defaults["energysim"] = strth.c_str(); 1342 query6.template_defaults["theta"] = strzd.c_str(); 1343 query6.template_defaults["spectrum"] = spectrum.c_str(); 1344 query6.template_defaults["estimator"] = estimator.c_str(); 1345 1346 if (print_queries) 1347 PrintQuery(query6.str()); 1348 1349 qlog << query6 << ";\n" << endl; 1350 if (connection.connected()) 1351 { 1352 cout << query6.execute().info() << endl; 1353 ShowWarnings(connection); 1354 Dump(flog, connection, "AnalysisMC"); 1355 1356 const auto sec6 = Time().UnixTime()-start6.UnixTime(); 1357 cout << "Execution time: " << sec6 << "s\n\n"; 1358 } 1359 1360 // ------------------------------------------------------------------- 1361 // ------------------------- SimulatedSpectrum ----------------------- 1362 // ------------------------------------------------------------------- 1363 1364 cout << separator("SimulatedSpectrum") << '\n'; 1365 1366 // [lo^(1+gammaa) - hi^(1+gamma)] / (1+gamma) 1367 1368 Time start7; 1369 /* 02:get-corsika-events.sql */ 1370 1371 // FIXME: Theta weights? 1372 // FIXME: energysim binning 1373 mysqlpp::Query query7(&connection); 1374 query7 << 1375 "CREATE TEMPORARY TABLE SimulatedSpectrum\n" 1376 "(\n" 1377 " `.energy` SMALLINT UNSIGNED NOT NULL COMMENT 'Bin Index [MC Energy]',\n" 1378 " CountN DOUBLE NOT NULL,\n" // COMMENT 'Number of events',\n" 1379 " CountW DOUBLE NOT NULL,\n" // COMMENT 'Sum of weights, reweighted sum',\n" 1380 " CountW2 DOUBLE NOT NULL,\n" // COMMENT 'Sum of square of weights'\n" 1381 " PRIMARY KEY (`.energy`) USING HASH\n" 1382 ") ENGINE=Memory\n" 1383 "AS\n" 1384 "(\n" 1385 " SELECT\n" 1386 " INTERVAL(LOG10(Energy), %100:energyest) AS `.energy`,\n" 1387 " COUNT(*) AS CountN,\n" 1388 " SUM( (%101:spectrum)/pow(Energy, SpectralIndex) ) AS CountW,\n" // Contents is: CountW 1389 " SUM(POW((%101:spectrum)/pow(Energy, SpectralIndex),2)) AS CountW2\n" // Error is: SQRT(CountW2) 1390 " FROM\n" 1391 " MonteCarloFiles\n" 1392 " LEFT JOIN\n" 1393 " factmc.RunInfoMC USING (FileId)\n" 1394 " LEFT JOIN\n" 1395 " factmc.OriginalMC USING (FileId)\n" 1396 " GROUP BY\n" 1397 " `.energy`\n" 1398 " ORDER BY\n" 1399 " `.energy`\n" 1400 ")"; 1401 1402 query7.parse(); 1403 for (auto it=env.cbegin(); it!=env.cend(); it++) 1404 query7.template_defaults[it->first.c_str()] = it->second.c_str(); 1405 1406 //query7.template_defaults["list"] = listmc.c_str(); 1407 query7.template_defaults["energyest"] = stre.c_str(); 1408 query7.template_defaults["spectrum"] = spectrum.c_str(); 1409 1410 if (print_queries) 1411 PrintQuery(query7.str()); 1412 1413 qlog << query7 << ";\n" << endl; 1414 if (connection.connected()) 1415 { 1416 cout << query7.execute().info() << endl; 1417 ShowWarnings(connection); 1418 Dump(flog, connection, "SimulatedSpectrum"); 1419 1420 const auto sec7 = Time().UnixTime()-start7.UnixTime(); 1421 cout << "Execution time: " << sec7 << "s\n"; 1422 1423 1424 if (verbose>0) 1425 { 1426 const mysqlpp::StoreQueryResult res7 = connection.query("SELECT * FROM SimulatedSpectrum").store(); 1427 1428 cout << " Bin CountW CountW2" << endl; 1429 const auto bins = binning_esim.vec(); 1430 for (auto ir=res7.cbegin(); ir!=res7.cend(); ir++) 1431 { 1432 const mysqlpp::Row &row = *ir; 1433 1434 const uint32_t bin = row[".energy"]; 1435 cout << setw(5) << bins[bin] << ": " << setw(10) << row["CountW"] << " " << setw(10) << row["CountW2"] << '\n'; 1436 } 1437 cout << endl; 1438 } 1439 } 1440 1441 1442 // ------------------------------------------------------------------- 1443 // ----------------------------- Spectrum ---------------------------- 1444 // ------------------------------------------------------------------- 1445 1446 cout << separator("Spectrum") << '\n'; 1447 1448 Time start8; 1449 1450 /* 02:calculate-spectrum.sql */ 1451 1452 mysqlpp::Query query8(&connection); 1453 query8 << 1454 // SELECT SUM(CountN) AS CountN, SUM(OnTime) AS OnTime FROM ThetaHist 1455 "CREATE TEMPORARY TABLE Spectrum\n" 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" 1456 1782 "(\n" 1457 1783 " `.energy` SMALLINT UNSIGNED NOT NULL COMMENT 'Bin Index [Energy]' PRIMARY KEY,\n" … … 1473 1799 " ErrSignalW DOUBLE NOT NULL COMMENT 'Error of weighted number of signal events in simulated data',\n" 1474 1800 " ErrBackgroundW DOUBLE NOT NULL COMMENT 'Error of weighted number of background events in simulated data',\n" 1475 " Flux DOUBLE NOT NULL COMMENT 'dN/dA/dt [cm^-2 s-^1 TeV^-1]',\n" 1476 " ErrFlux DOUBLE NOT NULL COMMENT 'dN/dA/dt [cm^-2 s-^1 TeV^-1]',\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" 1477 1807 " Bias DOUBLE NOT NULL COMMENT 'Energy Bias, average residual in lg(E)',\n" 1478 1808 " Resolution DOUBLE NOT NULL COMMENT 'Energy resolution, standard divation of residual in lg(E)',\n" 1479 " EfficiencyN DOUBLE NOT NULL COMMENT 'Simulated cut efficiency (weighted)',\n"1480 " EfficiencyW DOUBLE NOT NULL COMMENT 'Simulated cut efficiency (unweighted)',\n"1481 " ErrEfficiencyN DOUBLE NOT NULL COMMENT 'Error of simulated cut efficiency (weighted)',\n"1482 " ErrEfficiencyW DOUBLE NOT NULL COMMENT 'Error of simulated cut efficiency (unweighted)'\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" 1483 1813 ") ENGINE=Memory\n" 1484 "AS\n" 1485 "(\n" 1486 " WITH ThetaSums AS\n" 1487 " (\n" 1488 " SELECT\n" 1489 " SUM(CountN) AS CountSim,\n" 1490 " SUM(OnTime) AS ObsTime\n" 1491 " FROM\n" 1492 " ThetaDistribution\n" 1493 " ),\n" 1494 " Area AS\n" 1495 " (\n" 1496 " SELECT\n" 1497 " POW(MinImpactHi,2)*PI() AS Area\n" 1498 " FROM\n" 1499 " MonteCarloArea\n" 1500 " ),\n" 1501 " ResultMC AS\n" // Group AnalysisMC by EnergyEst Binning 1502 " (\n" 1503 " SELECT\n" 1504 " `.energyest` AS `.energy`,\n" 1505 " ANY_VALUE(SignalW) AS SignalW,\n" 1506 " ANY_VALUE(SignalW2) AS SignalW2,\n" 1507 " ANY_VALUE(BackgroundW) AS BackgroundW,\n" 1508 " ANY_VALUE(BackgroundW2) AS BackgroundW2,\n" 1509 " ANY_VALUE(SignalN) AS SignalN,\n" 1510 " ANY_VALUE(BackgroundN) AS BackgroundN,\n" 1511 " ANY_VALUE(ExcessW) AS ExcessW,\n" 1512 " ANY_VALUE(ExcessN) AS ExcessN,\n" 1513 " ANY_VALUE(ErrExcessW) AS ErrExcessW,\n" 1514 " ANY_VALUE(ErrExcessN) AS ErrExcessN,\n" 1515 " ANY_VALUE(BiasEst) AS Bias,\n" 1516 " ANY_VALUE(ResolutionEst) AS Resolution\n" 1517 " FROM\n" 1518 " AnalysisMC\n" 1519 " GROUP BY\n" 1520 " `.energy`\n" 1521 " ORDER BY\n" 1522 " `.energy`\n" 1523 " )\n" 1524 " SELECT\n" 1525 " `.energy`, lo, hi,\n" // Scale for Theta-Weights 1526 " `Signal`, `Background`/5 AS `Background`, `Excess`, `ErrExcess`, `Significance`,\n" 1527 " SQRT(`Signal`) AS ErrSignal,\n" 1528 " SQRT(`SignalW2`) AS ErrSignalW,\n" 1529 " SQRT(`Background`)/5 AS ErrBackground,\n" 1530 " SQRT(`BackgroundW2`)/5 AS ErrBackgroundW,\n" 1531 " ExcessN, ExcessW, ErrExcessN, ErrExcessW, SignalW, BackgroundW,\n" 1532 " AnalysisData.Excess/ResultMC.ExcessW*SimulatedSpectrum.CountW * 1000/(POW(10,hi)-POW(10,lo)) /Area/ObsTime / CountSim*ObsTime AS Flux,\n" 1533 " AnalysisData.Excess/ResultMC.ExcessW*SimulatedSpectrum.CountW * 1000/(POW(10,hi)-POW(10,lo)) /Area/ObsTime / CountSim*ObsTime\n" 1534 " * SQRT(\n" 1535 " + POW(AnalysisData.ErrExcess / AnalysisData.Excess, 2)\n" 1536 " + POW(ResultMC.ErrExcessW / ResultMC.ExcessW, 2)\n" 1537 " + SimulatedSpectrum.CountW2 / POW(SimulatedSpectrum.CountW,2)\n" 1538 " ) AS ErrFlux,\n" 1539 " Bias,\n" 1540 " Resolution,\n" 1541 " ResultMC.ExcessW/SimulatedSpectrum.CountW * CountSim/ObsTime AS EfficiencyW,\n" 1542 " ResultMC.ExcessN/SimulatedSpectrum.CountN AS EfficiencyN,\n" 1543 " ( POW(ResultMC.ErrExcessW/ResultMC.ExcessW, 2) + POW(SQRT(SimulatedSpectrum.CountW2)/SimulatedSpectrum.CountW, 2) )\n" 1544 " * POW(ResultMC.ExcessW/SimulatedSpectrum.CountW, 2) * CountSim/ObsTime AS ErrEfficiencyW,\n" 1545 " ( POW(ResultMC.ErrExcessN, 2) + POW(ResultMC.ExcessN, 2)/SimulatedSpectrum.CountN)/POW(SimulatedSpectrum.CountN, 2) AS ErrEfficiencyN,\n" 1546 " Area/10000*ResultMC.ExcessW/SimulatedSpectrum.CountW * CountSim/ObsTime AS EffAreaW,\n" 1547 " Area/10000*ResultMC.ExcessN/SimulatedSpectrum.CountN AS EffAreaN,\n" 1548 " Area/10000*( POW(ResultMC.ErrExcessW/ResultMC.ExcessW, 2) + POW(SQRT(SimulatedSpectrum.CountW2)/SimulatedSpectrum.CountW, 2) )\n" 1549 " * POW(ResultMC.ExcessW/SimulatedSpectrum.CountW, 2) * CountSim/ObsTime AS ErrEffAreaW,\n" 1550 " Area/10000*( POW(ResultMC.ErrExcessN, 2) + POW(ResultMC.ExcessN, 2)/SimulatedSpectrum.CountN)/POW(SimulatedSpectrum.CountN, 2) AS ErrEffAreaN\n" 1551 " FROM\n" 1552 " AnalysisData\n" 1553 " INNER JOIN\n" 1554 " ResultMC USING(`.energy`)\n" 1555 " INNER JOIN\n" 1556 " SimulatedSpectrum USING(`.energy`)\n" 1557 " INNER JOIN\n" 1558 " BinningEnergyEst ON `.energy`=bin\n" 1559 " CROSS JOIN\n" 1560 " ThetaSums, Area\n" 1561 " WHERE\n" 1562 " AnalysisData.Excess>0\n" 1563 " ORDER BY\n" 1564 " `.energy`\n" 1565 ")" 1566 ; 1567 1568 // [ Sa^2/a^2 + Sb^2/b^2 ] * a^2/b^2 1569 // [ (sc^2)/c^2+(sb^2)/b^2+(sa^2)/a^2 ] * a^2*b^2/c^2 1570 1571 query8.parse(); 1572 for (auto it=env.cbegin(); it!=env.cend(); it++) 1573 query8.template_defaults[it->first.c_str()] = it->second.c_str(); 1574 1575 //query8.template_defaults["area"] = area; 1576 //query8.template_defaults["ontime"] = resX[0]["OnTime"].data(); 1577 //query8.template_defaults["count"] = resX[0]["CountN"].data(); 1578 1579 if (print_queries) 1580 PrintQuery(query8.str()); 1581 1582 qlog << query8 << ";\n" << endl; 1583 if (connection.connected()) 1584 { 1585 cout << query8.execute().info() << endl; 1586 ShowWarnings(connection); 1587 Dump(flog, connection, "Spectrum"); 1588 1589 const auto sec8 = Time().UnixTime()-start8.UnixTime(); 1590 cout << "Execution time: " << sec8 << "s\n"; 1591 1592 1593 const mysqlpp::StoreQueryResult res8 = connection.query("SELECT * FROM Spectrum").store(); 1594 1595 if (verbose>0) 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") 1596 1838 { 1597 cout << " Bin Flux Error" << endl; 1598 const auto bins = binning_eest.vec(); 1599 for (auto ir=res8.cbegin(); ir!=res8.cend(); ir++) 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 1891 for (auto ir=res13.cbegin(); ir!=res13.cend(); ir++) 1600 1892 { 1893 // This is not efficient but easier to read and safer 1601 1894 const mysqlpp::Row &row = *ir; 1602 1895 1603 const uint32_t bin = row[".energy"]; 1604 cout << setw(5) << bins[bin] << ": " << setw(10) << row["Flux"] << " " << setw(10) << row["ErrFlux"] << '\n'; 1896 const size_t bin = row[*ib=="Theta" ? ".theta" : ".sparse_est"]; 1897 const double flux = row["Flux"]; 1898 const double error = row["ErrFlux"]; 1899 1900 #ifdef HAVE_ROOT 1901 const double dat_sig = row["Signal"]; 1902 const double dat_bg = row["Background"]; 1903 1904 const double eff = row["Efficiency"]; 1905 const double scale = row["Scale"]; 1906 1907 fc.SetMuMax(10*(dat_sig+dat_bg/5)); //has to be higher than sig+bg!! 1908 rolke.SetPoissonBkgKnownEff(dat_sig, dat_bg, 5, eff); 1909 1910 if (feldman) 1911 feldman_ul[bin] = fc.CalculateUpperLimit(dat_sig, dat_bg/5)*scale/eff; 1912 1913 rolke_ll[bin] = rolke.GetLowerLimit()*scale; 1914 rolke_ul[bin] = rolke.GetUpperLimit()*scale; 1915 #endif 1916 if (verbose>0) 1917 { 1918 cout << setw(5) << row["center"] << ": "; 1919 cout << setw(10) << row["Excess"] << " "; 1920 cout << setw(10) << flux << " "; 1921 cout << setw(10) << error << " "; 1922 cout << setw(10) << row["Significance"] << '\n'; 1923 } 1605 1924 } 1606 cout << endl; 1925 1926 Histogram hist_zd; 1927 hist_zd.dir = "Data/"+*ib; 1928 hist_zd.table = table; 1929 hist_zd.binningx = *ib=="Theta" ? binning_theta : binning_sparse; 1930 hist_zd.x = *ib=="Theta" ? ".theta" : ".sparse_est"; 1931 hist_zd.axisx = *ib=="Theta" ? "Zenith Distance #theta [#circ]" : "Energy lg(E/GeV)"; 1932 hist_zd.stats = false; 1933 1934 hist_zd.title = ""; 1935 hist_zd.v = "Signal"; 1936 hist_zd.err = "ErrSignal"; 1937 hist_zd.axisy = "Counts"; 1938 WriteHistogram(connection, hist_zd); 1939 1940 hist_zd.title = ""; 1941 hist_zd.v = "Background"; 1942 hist_zd.err = "ErrBackground"; 1943 hist_zd.axisy = "Counts"; 1944 WriteHistogram(connection, hist_zd); 1945 1946 hist_zd.title = ""; 1947 hist_zd.v = "Excess"; 1948 hist_zd.err = "ErrExcess"; 1949 hist_zd.axisy = "Counts"; 1950 WriteHistogram(connection, hist_zd); 1951 1952 hist_zd.title = ""; 1953 hist_zd.v = "Significance"; 1954 hist_zd.err = ""; 1955 hist_zd.axisy = "#sigma"; 1956 WriteHistogram(connection, hist_zd); 1957 1958 hist_zd.title = ""; 1959 hist_zd.v = "AvgEnergyEst"; 1960 hist_zd.err = ""; 1961 hist_zd.axisy = "<E_{est}>/GeV"; 1962 WriteHistogram(connection, hist_zd); 1963 1964 hist_zd.v = "ExcessRatio"; 1965 hist_zd.err = "ErrExcessRatio"; 1966 hist_zd.axisy = "Ratio"; 1967 WriteHistogram(connection, hist_zd); 1968 1969 1970 hist_zd.axisy = *ib=="Theta" ? "dN/dE [cm^{-2} s^{-1}]" : "dN/dE [cm^{-2} s^{-1} TeV^{-1}]"; 1971 1972 hist_zd.name = "Spectrum"; 1973 hist_zd.v = "Flux"; 1974 hist_zd.err = "ErrFlux"; 1975 WriteHistogram(connection, hist_zd); 1976 1977 #ifdef HAVE_ROOT 1978 hist_zd.axisy = *ib=="Theta" ? "UL [cm^{-2} s^{-1}]" : "UL [cm^{-2} s^{-1} TeV^{-1}]"; 1979 1980 if (feldman) 1981 { 1982 hist_zd.name = "FeldmanCousins"; 1983 WriteHistogram(connection, hist_zd, feldman_ul); 1984 } 1985 1986 hist_zd.name = "RolkeUL"; 1987 WriteHistogram(connection, hist_zd, rolke_ul); 1988 1989 hist_zd.name = "RolkeLL"; 1990 WriteHistogram(connection, hist_zd, rolke_ll); 1991 #endif 1607 1992 } 1608 1609 // -------------------------------------------------------------------------- 1610 1611 // Crab Nebula: 1 TeV: 3e-7 / m^2 / s / TeV 1612 // Crab Nebula: 1 TeV: 3e-11 / cm^2 / s / TeV 1613 1614 sindent indentm(mlog); 1615 1616 mlog << "void spectrum()\n"; 1617 mlog << "{\n" << indent(4); 1618 mlog << "// Energy Spectrum (e.g. Crab: 3e-11 [cm^-2 s-^1 TeV^-1] @ 1TeV)\n"; 1619 mlog << "TGraphErrors g;\n"; 1620 mlog << "g.SetNameTitle(\"Spectrum\", \"Energy Spectrum\");\n"; 1621 for (auto ir=res8.cbegin(); ir!=res8.cend(); ir++) 1622 { 1623 // This is not efficient but easier to read and safer 1624 const mysqlpp::Row &row = *ir; 1625 1626 const double hi = row["hi"]; 1627 const double lo = row["lo"]; 1628 const double center = (hi+lo)/2; 1629 const double flux = row["Flux"]; 1630 const double error = row["ErrFlux"]; 1631 1632 mlog << "g.SetPoint(g.GetN(), " << setw(8) << center << ", " << flux << ");\n"; 1633 mlog << "g.SetPointError(g.GetN()-1, 0, " << error << ");\n"; 1634 } 1635 mlog << "g.GetXaxis()->SetTitle(\"lg(E/TeV)\");\n"; 1636 mlog << "g.GetYaxis()->SetTitle(\"dN/dE [cm^{-2} s^{-1} TeV^{-1}]\");"; 1637 mlog << endl; 1638 } 1639 1640 WriteHistogram(connection, { 1641 .name = "Signal", 1642 .title = "Signal", 1643 .dir = "Eest/Measurement", 1644 .binningx = binning_eest, 1645 .binningy = {}, 1646 .table = "Spectrum", 1647 .x = ".energy", 1648 .y = "", 1649 .v = "Signal", 1650 .err = "ErrSignal", 1651 .axisx = "Energy lg(E_{est}/GeV)", 1652 .axisy = "Counts", 1653 .axisz = "", 1654 .stats = true 1655 }); 1656 1657 WriteHistogram(connection, { 1658 .name = "Background", 1659 .title = "Background", 1660 .dir = "Eest/Measurement", 1661 .binningx = binning_eest, 1662 .binningy = {}, 1663 .table = "Spectrum", 1664 .x = ".energy", 1665 .y = "", 1666 .v = "Background", 1667 .err = "ErrBackground", 1668 .axisx = "Energy lg(E_{est}/GeV)", 1669 .axisy = "Count", 1670 .axisz = "", 1671 .stats = true 1672 }); 1673 1674 WriteHistogram(connection, { 1675 .name = "Excess", 1676 .title = "Excess", 1677 .dir = "Eest/Measurement", 1678 .binningx = binning_eest, 1679 .binningy = {}, 1680 .table = "Spectrum", 1681 .x = ".energy", 1682 .y = "", 1683 .v = "Excess", 1684 .err = "ErrExcess", 1685 .axisx = "Energy lg(E_{est}/GeV)", 1686 .axisy = "Signal - Background (Counts)", 1687 .axisz = "", 1688 .stats = true 1689 }); 1690 1691 WriteHistogram(connection, { 1692 .name = "SignalW", 1693 .title = "SignalW", 1694 .dir = "Eest/Simulation/Weighted", 1695 .binningx = binning_eest, 1696 .binningy = {}, 1697 .table = "Spectrum", 1698 .x = ".energy", 1699 .y = "", 1700 .v = "SignalW", 1701 .err = "ErrSignalW", 1702 .axisx = "Energy lg(E_{est}/GeV)", 1703 .axisy = "Weighted", 1704 .axisz = "", 1705 .stats = true 1706 }); 1707 1708 WriteHistogram(connection, { 1709 .name = "BackgroundW", 1710 .title = "BackgroundW", 1711 .dir = "Eest/Simulation/Weighted", 1712 .binningx = binning_eest, 1713 .binningy = {}, 1714 .table = "Spectrum", 1715 .x = ".energy", 1716 .y = "", 1717 .v = "BackgroundW", 1718 .err = "ErrBackgroundW", 1719 .axisx = "Energy lg(E_{est}/GeV)", 1720 .axisy = "Weighted", 1721 .axisz = "", 1722 .stats = true 1723 }); 1724 1725 WriteHistogram(connection, { 1726 .name = "ExcessW", 1727 .title = "ExcessW", 1728 .dir = "Eest/Simulation/Weighted", 1729 .binningx = binning_eest, 1730 .binningy = {}, 1731 .table = "Spectrum", 1732 .x = ".energy", 1733 .y = "", 1734 .v = "ExcessW", 1735 .err = "ErrExcessW", 1736 .axisx = "Energy lg(E_{est}/GeV)", 1737 .axisy = "Signal - Background (Weighted)", 1738 .axisz = "", 1739 .stats = true 1740 }); 1741 1742 WriteHistogram(connection, { 1743 .name = "Significance", 1744 .title = "Significance", 1745 .dir = "Eest/Measurement", 1746 .binningx = binning_eest, 1747 .binningy = {}, 1748 .table = "Spectrum", 1749 .x = ".energy", 1750 .y = "", 1751 .v = "Significance", 1752 .err = "", 1753 .axisx = "Energy lg(E_{est}/GeV)", 1754 .axisy = "Li/Ma Significance", 1755 .axisz = "", 1756 .stats = true 1757 }); 1758 1759 WriteHistogram(connection, { 1760 .name = "Bias", 1761 .title = "Energy Bias", 1762 .dir = "Eest", 1763 .binningx = binning_eest, 1764 .binningy = {}, 1765 .table = "Spectrum", 1766 .x = ".energy", 1767 .y = "", 1768 .v = "Bias", 1769 .err = "Resolution", 1770 .axisx = "Energy lg(E_{sim}/GeV)", 1771 .axisy = "lg(E_{est}/E_{sim})", 1772 .axisz = "", 1773 .stats = false 1774 }); 1775 1776 WriteHistogram(connection, { 1777 .name = "Resolution", 1778 .title = "Energy Resolution", 1779 .dir = "Eest", 1780 .binningx = binning_eest, 1781 .binningy = {}, 1782 .table = "Spectrum", 1783 .x = ".energy", 1784 .y = "", 1785 .v = "Resolution", 1786 .err = "", 1787 .axisx = "Energy lg(E_{sim}/GeV)", 1788 .axisy = "#sigma(lg(E_{est}/E_{sim}))", 1789 .axisz = "", 1790 .stats = false 1791 }); 1792 1793 WriteHistogram(connection, { 1794 .name = "EfficiencyN", 1795 .title = "Efficiency (Counts)", 1796 .dir = "Eest/Efficiency", 1797 .binningx = binning_eest, 1798 .binningy = {}, 1799 .table = "Spectrum", 1800 .x = ".energy", 1801 .y = "", 1802 .v = "EfficiencyN", 1803 .err = "ErrEfficiencyN", 1804 .axisx = "Energy lg(E_{sim}/GeV)", 1805 .axisy = "Efficiency", 1806 .axisz = "", 1807 .stats = true 1808 }); 1809 1810 WriteHistogram(connection, { 1811 .name = "EfficiencyW", 1812 .title = "Efficiency (Weighted)", 1813 .dir = "Eest/Efficiency", 1814 .binningx = binning_eest, 1815 .binningy = {}, 1816 .table = "Spectrum", 1817 .x = ".energy", 1818 .y = "", 1819 .v = "EfficiencyW", 1820 .err = "ErrEfficiencyW", 1821 .axisx = "Energy lg(E_{sim}/GeV)", 1822 .axisy = "Efficiency", 1823 .axisz = "", 1824 .stats = true 1825 }); 1826 1827 WriteHistogram(connection, { 1828 .name = "EffAreaN", 1829 .title = "Effective Area (Counts)", 1830 .dir = "Eest/EffArea", 1831 .binningx = binning_eest, 1832 .binningy = {}, 1833 .table = "Spectrum", 1834 .x = ".energy", 1835 .y = "", 1836 .v = "EffAreaN", 1837 .err = "ErrEffAreaN", 1838 .axisx = "Energy lg(E_{sim}/GeV)", 1839 .axisy = "Effective Area A [m^{2}]", 1840 .axisz = "", 1841 .stats = false 1842 }); 1843 1844 WriteHistogram(connection, { 1845 .name = "EffAreaW", 1846 .title = "Effective Area (Weighted)", 1847 .dir = "Eest/EffArea", 1848 .binningx = binning_eest, 1849 .binningy = {}, 1850 .table = "Spectrum", 1851 .x = ".energy", 1852 .y = "", 1853 .v = "EffAreaW", 1854 .err = "ErrEffAreaW", 1855 .axisx = "Energy lg(E_{sim}/GeV)", 1856 .axisy = "Effective Area A [m^{2}]", 1857 .axisz = "", 1858 .stats = false 1859 }); 1860 1861 WriteHistogram(connection, { 1862 .name = "Spectrum", 1863 .title = "Differential Energy Spectrum", 1864 .dir = "", 1865 .binningx = binning_eest, 1866 .binningy = {}, 1867 .table = "Spectrum", 1868 .x = ".energy", 1869 .y = "", 1870 .v = "Flux", 1871 .err = "ErrFlux", 1872 .axisx = "Energy lg(E/GeV)", 1873 .axisy = "dN/dE [cm^{-2} s^{-1} TeV^{-1}]", 1874 .axisz = "", 1875 .stats = false 1876 }); 1877 1878 1879 // ------------------------------------------------------------------- 1880 // ---------------------------- Threshold ---------------------------- 1881 // ------------------------------------------------------------------- 1882 1883 cout << separator("Threshold") << '\n'; 1884 1885 Time start9; 1886 1887 /* 02:calculate-threshold.sql */ 1888 1889 // This query gets the analysis results versus Simulated Energy from the combined table 1890 mysqlpp::Query query9(&connection); 1891 query9 << 1892 // SELECT SUM(CountN) AS CountN, SUM(OnTime) AS OnTime FROM ThetaHist 1893 "CREATE TEMPORARY TABLE Threshold ENGINE=Memory AS\n" 1894 "(\n" 1895 " WITH ThetaSums AS\n" 1896 " (\n" 1897 " SELECT\n" 1898 " SUM(CountN) AS CountSim,\n" 1899 " SUM(OnTime) AS ObsTime\n" 1900 " FROM\n" 1901 " ThetaDistribution\n" 1902 " ),\n" 1903 " Area AS\n" 1904 " (\n" 1905 " SELECT\n" 1906 " POW(MinImpactHi,2)*PI() AS Area\n" 1907 " FROM\n" 1908 " MonteCarloArea\n" 1909 " ),\n" 1910 " ResultMC AS\n" // Group AnalysisMC by EnergySim Binning 1911 " (\n" 1912 " SELECT\n" 1913 " `.energysim` AS `.energy`,\n" 1914 " ANY_VALUE(ThresholdW) AS ThresholdW,\n" 1915 " ANY_VALUE(ThresholdW2) AS ThresholdW2,\n" 1916 " ANY_VALUE(ThresholdN) AS ThresholdN,\n" 1917 " ANY_VALUE(BiasSim) AS Bias,\n" 1918 " ANY_VALUE(ResolutionSim) AS Resolution\n" 1919 " FROM\n" 1920 " AnalysisMC\n" 1921 " GROUP BY\n" 1922 " `.energy`\n" 1923 " )\n" 1924 " SELECT\n" 1925 " `.energy`, lo, hi,\n" 1926 " ThresholdW,\n" 1927 " SQRT(ThresholdW2) AS ErrThresholdW,\n" 1928 " ThresholdN,\n" 1929 " SQRT(ThresholdN) AS ErrThresholdN,\n" // Scale for Theta-Weights 1930 " ThresholdW * 1000/(POW(10,hi)-POW(10,lo)) / Area / CountSim*ObsTime AS Flux,\n" 1931 " SQRT(ThresholdW2) * 1000/(POW(10,hi)-POW(10,lo)) / Area / CountSim*ObsTime AS ErrFlux,\n" 1932 " ThresholdN/CountN AS CutEfficiencyN,\n" 1933 " ThresholdW/CountW AS CutEfficiencyW,\n" 1934 " ThresholdN/CountN*(1/ThresholdN + 1/CountN) AS ErrCutEfficiencyN,\n" 1935 " ThresholdW/CountW*(ThresholdW2/POW(ThresholdW,2) + CountW2/POW(CountW,2)) AS ErrCutEfficiencyW,\n" 1936 " Bias,\n" 1937 " Resolution\n" 1938 " FROM\n" 1939 " ResultMC\n" 1940 " INNER JOIN\n" 1941 " Triggered USING (`.energy`)\n" 1942 " INNER JOIN\n" 1943 " BinningEnergySim ON `.energy`=bin\n" 1944 " CROSS JOIN\n" 1945 " ThetaSums, Area\n" 1946 " WHERE\n" 1947 " ThresholdW>0 AND ThresholdW2>0\n" 1948 " ORDER BY\n" 1949 " `.energy`\n" 1950 ")"; 1951 1952 query9.parse(); 1953 for (auto it=env.cbegin(); it!=env.cend(); it++) 1954 query9.template_defaults[it->first.c_str()] = it->second.c_str(); 1955 1956 //query9.template_defaults["area"] = area; 1957 //query9.template_defaults["ontime"] = resX[0]["OnTime"].data(); 1958 //query9.template_defaults["count"] = resX[0]["CountN"].data(); 1959 1960 if (print_queries) 1961 PrintQuery(query9.str()); 1962 1963 qlog << query9 << ";\n" << endl; 1964 if (connection.connected()) 1965 { 1966 cout << query9.execute().info() << endl; 1967 ShowWarnings(connection); 1968 Dump(flog, connection, "Threshold"); 1969 1970 const auto sec9 = Time().UnixTime()-start9.UnixTime(); 1971 cout << "Execution time: " << sec9 << "s\n"; 1972 1973 // -------------------------------------------------------------------------- 1974 1975 const mysqlpp::StoreQueryResult res9 = connection.query("SELECT * FROM Threshold").store(); 1976 1977 sindent indentm(mlog, 4); 1978 1979 mlog << '\n'; 1980 mlog << "// Energy Threshold\n"; 1981 mlog << "TGraphErrors g;\n"; 1982 mlog << "g.SetNameTitle(\"Threshold\", \"Energy Threshold\");\n"; 1983 for (auto ir=res9.cbegin(); ir!=res9.cend(); ir++) 1984 { 1985 // This is not efficient but easier to read and safer 1986 const mysqlpp::Row &row = *ir; 1987 1988 const double hi = row["hi"]; 1989 const double lo = row["lo"]; 1990 const double center = (hi+lo)/2; 1991 const double width = pow(10, hi)-pow(10, lo); 1992 const double flux = row["Flux"] /width; 1993 const double error = row["ErrFlux"]/width; 1994 1995 mlog << "g.SetPoint(g.GetN(), " << setw(8) << center << ", " << flux << ");\n"; 1996 mlog << "g.SetPointError(g.GetN()-1, 0, " << error << ");\n"; 1997 } 1998 mlog << "g.GetXaxis()->SetTitle(\"lg(E/TeV)\");\n"; 1999 mlog << "g.GetYaxis()->SetTitle(\"dN/dE [cm^{-2} s^{-1} TeV^{-1}]\");\n"; 2000 mlog << indent(0) << "}" << endl; 2001 } 2002 2003 WriteHistogram(connection, { 2004 .name = "SimExcessW", 2005 .title = "Weighted Simulated Excess", 2006 .dir = "Esim/Simulation/Weighted", 2007 .binningx = binning_esim, 2008 .binningy = {}, 2009 .table = "Threshold", 2010 .x = ".energy", 2011 .y = "", 2012 .v = "ThresholdW", 2013 .err = "ErrThresholdW", 2014 .axisx = "Energy lg(E_{sim}/GeV)", 2015 .axisy = "Weighted Counts", 2016 .axisz = "", 2017 .stats = true 2018 }); 2019 2020 WriteHistogram(connection, { 2021 .name = "SimExcessN", 2022 .title = "Simulated Excess", 2023 .dir = "Esim/Simulation/Counts", 2024 .binningx = binning_esim, 2025 .binningy = {}, 2026 .table = "Threshold", 2027 .x = ".energy", 2028 .y = "", 2029 .v = "ThresholdN", 2030 .err = "ErrThresholdN", 2031 .axisx = "Energy lg(E_{sim}/GeV)", 2032 .axisy = "Counts", 2033 .axisz = "", 2034 .stats = true 2035 }); 2036 2037 WriteHistogram(connection, { 2038 .name = "SimSpectrumW", 2039 .title = "Weighted Simulated Excess Spectrum", 2040 .dir = "Esim/Simulation/Weighted", 2041 .binningx = binning_esim, 2042 .binningy = {}, 2043 .table = "Threshold", 2044 .x = ".energy", 2045 .y = "", 2046 .v = "Flux", 2047 .err = "ErrFlux", 2048 .axisx = "Energy lg(E_{sim}/GeV)", 2049 .axisy = "dN/dE [cm^{-2} s^{-1} TeV^{-1}]", 2050 .axisz = "", 2051 .stats = true 2052 }); 2053 2054 WriteHistogram(connection, { 2055 .name = "CutEfficiencyN", 2056 .title = "Cut Efficency (Unweighted)", 2057 .dir = "Esim/CutEfficiency", 2058 .binningx = binning_esim, 2059 .binningy = {}, 2060 .table = "Threshold", 2061 .x = ".energy", 2062 .y = "", 2063 .v = "CutEfficiencyN", 2064 .err = "ErrCutEfficiencyN", 2065 .axisx = "Energy lg(E_{sim}/GeV)", 2066 .axisy = "Efficiency", 2067 .axisz = "", 2068 .stats = false 2069 }); 2070 2071 WriteHistogram(connection, { 2072 .name = "CutEfficiencyW", 2073 .title = "Cut Efficency (Weighted)", 2074 .dir = "Esim/CutEfficiency", 2075 .binningx = binning_esim, 2076 .binningy = {}, 2077 .table = "Threshold", 2078 .x = ".energy", 2079 .y = "", 2080 .v = "CutEfficiencyW", 2081 .err = "ErrCutEfficiencyW", 2082 .axisx = "Energy lg(E_{sim}/GeV)", 2083 .axisy = "Efficiency", 2084 .axisz = "", 2085 .stats = false 2086 }); 2087 2088 WriteHistogram(connection, { 2089 .name = "Bias", 2090 .title = "Energy Bias", 2091 .dir = "Esim", 2092 .binningx = binning_esim, 2093 .binningy = {}, 2094 .table = "Threshold", 2095 .x = ".energy", 2096 .y = "", 2097 .v = "Bias", 2098 .err = "Resolution", 2099 .axisx = "Energy lg(E_{sim}/GeV)", 2100 .axisy = "lg(E_{est}/E_{sim})", 2101 .axisz = "", 2102 .stats = false 2103 }); 2104 2105 WriteHistogram(connection, { 2106 .name = "Resolution", 2107 .title = "Energy Resolution", 2108 .dir = "Esim", 2109 .binningx = binning_esim, 2110 .binningy = {}, 2111 .table = "Threshold", 2112 .x = ".energy", 2113 .y = "", 2114 .v = "Resolution", 2115 .err = "", 2116 .axisx = "Energy lg(E_{sim}/GeV)", 2117 .axisy = "#sigma(lg(E_{est}/E_{sim}))", 2118 .axisz = "", 2119 .stats = false 2120 }); 2121 2122 2123 // ------------------------------------------------------------------- 2124 // ---------------------------- Migration ---------------------------- 2125 // ------------------------------------------------------------------- 2126 2127 cout << separator("Migration") << '\n'; 2128 2129 Time start10; 2130 2131 /* 02:obtain-migration-matrix.sql */ 2132 2133 // This query gets the analysis results versus Estimated Energy from the combined table 2134 mysqlpp::Query query10(&connection); 2135 query10 << 2136 "CREATE TEMPORARY TABLE Migration ENGINE=Memory AS\n" 2137 "(\n" 2138 " SELECT\n" 2139 " `.energyest`,\n" 2140 " `.energysim`,\n" 2141 " BinningEnergySim.lo AS EsimLo,\n" 2142 " BinningEnergySim.hi AS EsimHi,\n" 2143 " BinningEnergyEst.lo AS EestLo,\n" 2144 " BinningEnergyEst.hi AS EestHi,\n" 2145 " ANY_VALUE(MigrationW) AS MigrationW,\n" 2146 " ANY_VALUE(MigrationN) AS MigrationN\n" 2147 // FIXME: Errors 2148 " FROM\n" 2149 " AnalysisMC\n" 2150 " INNER JOIN\n" 2151 " BinningEnergyEst ON `.energyest`=BinningEnergyEst.bin\n" 2152 " INNER JOIN\n" 2153 " BinningEnergySim ON `.energysim`=BinningEnergySim.bin\n" 2154 " GROUP BY\n" 2155 " `.energyest`, `.energysim`\n" 2156 " ORDER BY\n" 2157 " `.energyest`, `.energysim`\n" 2158 ")"; 2159 2160 if (print_queries) 2161 PrintQuery(query10.str()); 2162 2163 qlog << query10 << ";\n" << endl; 2164 if (connection.connected()) 2165 { 2166 cout << query10.execute().info() << endl; 2167 ShowWarnings(connection); 2168 Dump(flog, connection, "Migration"); 2169 2170 const auto sec10 = Time().UnixTime()-start10.UnixTime(); 2171 cout << "Execution time: " << sec10 << "s\n\n"; 2172 } 2173 2174 WriteHistogram(connection, { 2175 .name = "MigrationN", 2176 .title = "Energy Migration", 2177 .dir = "", 2178 .binningx = binning_esim, 2179 .binningy = binning_eest, 2180 .table = "Migration", 2181 .x = ".energysim", 2182 .y = ".energyest", 2183 .v = "MigrationN", 2184 .err = "", 2185 .axisx = "Energy lg(E_{sim}/GeV)", 2186 .axisy = "Energy lg(E_{est}/GeV)", 2187 .axisz = "Counts", 2188 .stats = false 2189 }); 2190 2191 WriteHistogram(connection, { 2192 .name = "MigrationW", 2193 .title = "Energy Migration", 2194 .dir = "", 2195 .binningx = binning_esim, 2196 .binningy = binning_eest, 2197 .table = "Migration", 2198 .x = ".energysim", 2199 .y = ".energyest", 2200 .v = "MigrationW", 2201 .err = "", 2202 .axisx = "Energy lg(E_{sim}/GeV)", 2203 .axisy = "Energy lg(E_{est}/GeV)", 2204 .axisz = "Weighted Counts", 2205 .stats = false 2206 }); 2207 2208 2209 2210 // ------------------------------------------------------------------- 2211 // --------------------------- 11: Summary --------------------------- 1993 } 1994 1995 // ------------------------------------------------------------------- 1996 // ----------------------------- Summary ----------------------------- 1997 // ------------------------------------------------------------------- 2212 1998 2213 1999 cout << separator("Summary") << '\n';
Note:
See TracChangeset
for help on using the changeset viewer.
