Changeset 19889
 Timestamp:
 Dec 12, 2019, 4:25:49 PM (6 weeks ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/FACT++/src/spectrum.cc
r19888 r19889 691 691 692 692 //  693 //  0: Binnings  693 //  Binnings  694 //  694 695 695 696 cout << separator("Binnings") << '\n'; … … 704 705 705 706 //  706 //  1: DataFiles  707 //  DataFiles  708 //  709 707 710 cout << separator("DataFiles") << '\n'; 708 711 … … 754 757 755 758 //  756 //  2: ObservationTime  759 //  ObservationTime  760 //  757 761 758 762 cout << separator("ObservationTime") << '\n'; … … 867 871 868 872 //  869 //  3: MonteCarloFiles  873 //  MonteCarloFiles  874 //  870 875 871 876 cout << separator("MonteCarloFiles") << '\n'; … … 931 936 932 937 //  933 //  3b: Monte Carlo Area  938 //  Monte Carlo Area  939 //  934 940 935 941 cout << separator("MonteCarloArea") << '\n'; … … 982 988 983 989 //  984 //  4: ThetaHist  985 986 cout << separator("ThetaHist") << '\n'; 990 //  ThetaDistribution  991 //  992 993 cout << separator("ThetaDistribution") << '\n'; 987 994 988 995 Time start4; … … 991 998 mysqlpp::Query query4(&connection); 992 999 query4 << 993 "CREATE TEMPORARY TABLE Theta Hist\n"1000 "CREATE TEMPORARY TABLE ThetaDistribution\n" 994 1001 "(\n" 995 1002 " `.theta` SMALLINT UNSIGNED NOT NULL,\n" … … 1048 1055 cout << query4.execute().info() << endl; 1049 1056 ShowWarnings(connection); 1050 Dump(flog, connection, "Theta Hist");1057 Dump(flog, connection, "ThetaDistribution"); 1051 1058 1052 1059 const auto sec4 = Time().UnixTime()start4.UnixTime(); … … 1056 1063 if (verbose>0) 1057 1064 { 1058 const mysqlpp::StoreQueryResult res4 = connection.query("SELECT * FROM Theta Hist").store();1065 const mysqlpp::StoreQueryResult res4 = connection.query("SELECT * FROM ThetaDistribution").store(); 1059 1066 1060 1067 cout << " Bin MC Counts OnTime ZdWeight\n"; … … 1077 1084 .binningx = binning_theta, 1078 1085 .binningy = {}, 1079 .table = "Theta Hist",1086 .table = "ThetaDistribution", 1080 1087 .x = ".theta", 1081 1088 .y = "", … … 1094 1101 .binningx = binning_theta, 1095 1102 .binningy = {}, 1096 .table = "Theta Hist",1103 .table = "ThetaDistribution", 1097 1104 .x = ".theta", 1098 1105 .y = "", … … 1111 1118 .binningx = binning_theta, 1112 1119 .binningy = {}, 1113 .table = "Theta Hist",1120 .table = "ThetaDistribution", 1114 1121 .x = ".theta", 1115 1122 .y = "", … … 1125 1132 //  1126 1133 //  5: AnalysisData  1134 //  1127 1135 1128 1136 cout << separator("AnalysisData") << '\n'; … … 1201 1209 1202 1210 //  1203 //  6: ResultMC  1211 //  Triggered  1212 //  1213 1214 cout << separator("Triggered") << '\n'; 1215 1216 Time start6a; 1217 1218 /* 02:analyzesimulation.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 //  1204 1284 1205 1285 cout << separator("ResultMC") << '\n'; … … 1278 1358 } 1279 1359 1280 1281 //  1282 //  7: SimulatedSpectrum1360 //  1361 //  SimulatedSpectrum  1362 //  1283 1363 1284 1364 cout << separator("SimulatedSpectrum") << '\n'; … … 1361 1441 1362 1442 //  1363 //  8: Spectrum  1443 //  Spectrum  1444 //  1364 1445 1365 1446 cout << separator("Spectrum") << '\n'; … … 1409 1490 " SUM(OnTime) AS ObsTime\n" 1410 1491 " FROM\n" 1411 " Theta Hist\n"1492 " ThetaDistribution\n" 1412 1493 " ),\n" 1413 1494 " Area AS\n" … … 1462 1543 " ( POW(ResultMC.ErrExcessW/ResultMC.ExcessW, 2) + POW(SQRT(SimulatedSpectrum.CountW2)/SimulatedSpectrum.CountW, 2) )\n" 1463 1544 " * POW(ResultMC.ExcessW/SimulatedSpectrum.CountW, 2) * CountSim/ObsTime AS ErrEfficiencyW,\n" 1464 " ( POW(ResultMC.ErrExcessN, 2) + POW(ResultMC.ExcessN, 2)/SimulatedSpectrum.CountN)/POW(SimulatedSpectrum.CountN, 2) AS ErrEfficiencyN\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" 1465 1551 " FROM\n" 1466 1552 " AnalysisData\n" … … 1740 1826 1741 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, { 1742 1862 .name = "Spectrum", 1743 1863 .title = "Differential Energy Spectrum", … … 1758 1878 1759 1879 //  1760 //  9: Threshold  1880 //  Threshold  1881 //  1761 1882 1762 1883 cout << separator("Threshold") << '\n'; … … 1778 1899 " SUM(OnTime) AS ObsTime\n" 1779 1900 " FROM\n" 1780 " Theta Hist\n"1901 " ThetaDistribution\n" 1781 1902 " ),\n" 1782 1903 " Area AS\n" … … 1809 1930 " ThresholdW * 1000/(POW(10,hi)POW(10,lo)) / Area / CountSim*ObsTime AS Flux,\n" 1810 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" 1811 1936 " Bias,\n" 1812 1937 " Resolution\n" 1813 1938 " FROM\n" 1814 1939 " ResultMC\n" 1940 " INNER JOIN\n" 1941 " Triggered USING (`.energy`)\n" 1815 1942 " INNER JOIN\n" 1816 1943 " BinningEnergySim ON `.energy`=bin\n" … … 1926 2053 1927 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, { 1928 2089 .name = "Bias", 1929 2090 .title = "Energy Bias", … … 1961 2122 1962 2123 //  1963 //  10: Migration  2124 //  Migration  2125 //  1964 2126 1965 2127 cout << separator("Migration") << '\n';
Note: See TracChangeset
for help on using the changeset viewer.