Changeset 19898 for trunk/FACT++

Timestamp:

12/14/19 23:42:10 (5 years ago)

Author:

tbretz

Message:

Updated to include Integral and Integrated Results.

File:

• trunk/FACT++/src/spectrum.cc (modified) (6 diffs)

trunk/FACT++/src/spectrum.cc

@@ -1376 +1376 @@
         query9.template_defaults["binning"]  = *ib=="theta" ? "BinningTheta" : ("BinningEnergy_"+*ib).c_str();
         query9.template_defaults["bin"]      = *ib=="theta" ? "`.theta`"     : ("`."+*ib+"_sim`").c_str();
-        query9.template_defaults["binwidth"] = *ib=="theta" ? "1"            : "1000/(POW(10,hi)-POW(10,lo))";
+        query9.template_defaults["binwidth"] = *ib=="theta" ? "1"            : "(POW(10,hi)-POW(10,lo))/1000";
 
         if (print_queries)
@@ -1822 +1822 @@
             << indent(3) << spectrum_sql << indent(0) <<
             ")";
-*/
+
         // [ Sa^2/a^2 + Sb^2/b^2 ] * a^2/b^2
         // [ (sc^2)/c^2+(sb^2)/b^2+(sa^2)/a^2 ]  * a^2*b^2/c^2
@@ -1831 +1831 @@
 
         query13.template_defaults["table"]     = table.c_str();
-        query13.template_defaults["bin"]       = *ib=="Theta"       ? "`.theta`"    : "`.sparse_est`";
-        query13.template_defaults["id"]        = *ib=="Theta"       ? "Sim"         : "Est";
-        query13.template_defaults["weight"]    = *ib=="Theta"       ? "ZdWeight"    : "1";
-        query13.template_defaults["errweight"] = *ib=="Theta"       ? "ErrZdWeight" : "1";
+        query13.template_defaults["bin"]       = *ib=="Theta" ? "`.theta`"    : "`.sparse_est`";
+        query13.template_defaults["id"]        = *ib=="Theta" ? "Sim"         : "Est";
+        query13.template_defaults["weight"]    = *ib=="Theta" ? "ZdWeight"    : "1";
+        query13.template_defaults["errweight"] = *ib=="Theta" ? "ErrZdWeight" : "1";
         if (*ib=="Theta")
         {
@@ -1888 +1888 @@
             map<size_t, double> rolke_ul;
             map<size_t, double> rolke_ll;
+            map<size_t, double> rolke_int;
 
             for (auto ir=res13.cbegin(); ir!=res13.cend(); ir++)
@@ -1894 +1895 @@
                 const mysqlpp::Row &row = *ir;
 
-                const size_t bin     = row[*ib=="Theta" ? ".theta" : ".sparse_est"];
-                const double flux    = row["Flux"];
-                const double error   = row["ErrFlux"];
+                const size_t bin      = row[*ib=="Theta" ? ".theta" : ".sparse_est"];
+                const double flux     = row["Flux"];
+                const double error    = row["ErrFlux"];
 
 #ifdef HAVE_ROOT
-                const double dat_sig = row["Signal"];
-                const double dat_bg  = row["Background"];
-
-                const double eff     = row["Efficiency"];
-                const double scale   = row["Scale"];
+                const double dat_sig  = row["Signal"];
+                const double dat_bg   = row["Background"];
+
+                const double dat_isig = row["IntegralSignal"];
+                const double dat_ibg  = row["IntegralBackground"];
+
+                const double eff      = row["Efficiency"];
+                const double ieff     = row["IntegralEfficiency"];
+
+                const double areatime = row["AreaTime"];
+                const double width    = row["Width"];
 
                 fc.SetMuMax(10*(dat_sig+dat_bg/5)); //has to be higher than sig+bg!!
-                rolke.SetPoissonBkgKnownEff(dat_sig, dat_bg, 5, eff);
 
                 if (feldman)
-                    feldman_ul[bin] = fc.CalculateUpperLimit(dat_sig, dat_bg/5)*scale/eff;
-
-                rolke_ll[bin] = rolke.GetLowerLimit()*scale;
-                rolke_ul[bin] = rolke.GetUpperLimit()*scale;
+                    feldman_ul[bin] = fc.CalculateUpperLimit(dat_sig, dat_bg/5)/width/areatime/eff;
+
+                rolke.SetPoissonBkgKnownEff(dat_sig,  dat_bg, 5, eff);
+                rolke_ll[bin]  = rolke.GetLowerLimit()/width/areatime;
+                rolke_ul[bin]  = rolke.GetUpperLimit()/width/areatime;
+
+                rolke.SetPoissonBkgKnownEff(dat_isig, dat_ibg, 5, ieff);
+                rolke_int[bin] = rolke.GetUpperLimit()/areatime;
 #endif
                 if (verbose>0)
@@ -1924 +1934 @@
             }
 
-            Histogram hist_zd;
-            hist_zd.dir      = "Data/"+*ib;
-            hist_zd.table    = table;
-            hist_zd.binningx = *ib=="Theta" ? binning_theta : binning_sparse;
-            hist_zd.x        = *ib=="Theta" ? ".theta" : ".sparse_est";
-            hist_zd.axisx    = *ib=="Theta" ? "Zenith Distance #theta [#circ]" : "Energy lg(E/GeV)";
-            hist_zd.stats    = false;
-
-            hist_zd.title    = "";
-            hist_zd.v        = "Signal";
-            hist_zd.err      = "ErrSignal";
-            hist_zd.axisy    = "Counts";
-            WriteHistogram(connection, hist_zd);
-
-            hist_zd.title    = "";
-            hist_zd.v        = "Background";
-            hist_zd.err      = "ErrBackground";
-            hist_zd.axisy    = "Counts";
-            WriteHistogram(connection, hist_zd);
-
-            hist_zd.title    = "";
-            hist_zd.v        = "Excess";
-            hist_zd.err      = "ErrExcess";
-            hist_zd.axisy    = "Counts";
-            WriteHistogram(connection, hist_zd);
-
-            hist_zd.title    = "";
-            hist_zd.v        = "Significance";
-            hist_zd.err      = "";
-            hist_zd.axisy    = "#sigma";
-            WriteHistogram(connection, hist_zd);
-
-            hist_zd.title    = "";
-            hist_zd.v        = "AvgEnergyEst";
-            hist_zd.err      = "";
-            hist_zd.axisy    = "<E_{est}>/GeV";
-            WriteHistogram(connection, hist_zd);
-
-            hist_zd.v        = "ExcessRatio";
-            hist_zd.err      = "ErrExcessRatio";
-            hist_zd.axisy    = "Ratio";
-            WriteHistogram(connection, hist_zd);
-
-
-            hist_zd.axisy    = *ib=="Theta" ? "dN/dE [cm^{-2} s^{-1}]" : "dN/dE [cm^{-2} s^{-1} TeV^{-1}]";
-
-            hist_zd.name     = "Spectrum";
-            hist_zd.v        = "Flux";
-            hist_zd.err      = "ErrFlux";
-            WriteHistogram(connection, hist_zd);
+            Histogram hist;
+            hist.dir      = "Data/"+*ib;
+            hist.table    = table;
+            hist.binningx = *ib=="Theta" ? binning_theta : binning_sparse;
+            hist.x        = *ib=="Theta" ? ".theta" : ".sparse_est";
+            hist.axisx    = *ib=="Theta" ? "Zenith Distance #theta [#circ]" : "Energy lg(E/GeV)";
+            hist.stats    = false;
+
+            const vector<string> types { "", "Integral" };
+            for (auto it=types.cbegin(); it!=types.cend(); it++)
+            {
+                hist.axisy = "Counts";
+                if (*ib=="Energy")
+                    hist.axisy += " (E>E_{lo})";
+
+                hist.title = "";
+                hist.v     = *it+"Signal";
+                hist.err   = "Err"+*it+"Signal";
+                WriteHistogram(connection, hist);
+
+                hist.title = "";
+                hist.v     = *it+"Background";
+                hist.err   = "Err"+*it+"Background";
+                WriteHistogram(connection, hist);
+
+                hist.title = "";
+                hist.v     = *it+"Excess";
+                hist.err   = "Err"+*it+"Excess";
+                WriteHistogram(connection, hist);
+
+                hist.title = "";
+                hist.v     = *it+"Significance";
+                hist.err   = "";
+                hist.axisy = "#sigma";
+                if (*ib=="Energy")
+                    hist.axisy += " (E>E_{lo})";
+                WriteHistogram(connection, hist);
+
+                hist.title = "";
+                hist.v     = *it+"AvgEnergyEst";
+                hist.err   = "";
+                hist.axisy = "<E_{est}>/GeV";
+                if (*ib=="Energy")
+                    hist.axisy += " (E>E_{lo})";
+                WriteHistogram(connection, hist);
+
+                hist.title = "";
+                hist.v     = *it+"ExcessRatio";
+                hist.err   = "Err"+*it+"ExcessRatio";
+                hist.axisy = "Ratio";
+                if (*ib=="Energy")
+                    hist.axisy += " (E>E_{lo})";
+                WriteHistogram(connection, hist);
+            }
+
+            hist.axisy = "dN/dE ";
+            if (*ib=="Energy")
+                hist.axisy += "(E>E_{lo}) ";
+            hist.axisy += *ib=="Theta" ? "[cm^{-2} s^{-1}]" : "[cm^{-2} s^{-1} TeV^{-1}]";
+
+            hist.name = "Spectrum";
+            hist.v    = "Flux";
+            hist.err  = "ErrFlux";
+            WriteHistogram(connection, hist);
+
+            if (*ib=="Energy")
+            {
+                hist.axisy = "dN/dE (E>E_{lo}) [cm^{-2} s^{-1}]";
+
+                hist.name  = "IntegralSpectrum";
+                hist.v     = "IntegralFlux";
+                hist.err   = "ErrIntegralFlux";
+                WriteHistogram(connection, hist);
+
+                hist.name  = "IntegratedSpectrum";
+                hist.v     = "IntegratedFlux";
+                hist.err   = "ErrIntegratedFlux";
+                WriteHistogram(connection, hist);
+            }
 
 #ifdef HAVE_ROOT
-            hist_zd.axisy    = *ib=="Theta" ? "UL [cm^{-2} s^{-1}]" : "UL [cm^{-2} s^{-1} TeV^{-1}]";
+            hist.axisy = *ib=="Theta" ? "UL [cm^{-2} s^{-1}]" : "UL [cm^{-2} s^{-1} TeV^{-1}]";
 
             if (feldman)
             {
-                hist_zd.name = "FeldmanCousins";
-                WriteHistogram(connection, hist_zd, feldman_ul);
+                hist.name = "FeldmanCousins";
+                WriteHistogram(connection, hist, feldman_ul);
             }
 
-            hist_zd.name = "RolkeUL";
-            WriteHistogram(connection, hist_zd, rolke_ul);
-
-            hist_zd.name = "RolkeLL";
-            WriteHistogram(connection, hist_zd, rolke_ll);
+            hist.name = "RolkeUL";
+            WriteHistogram(connection, hist, rolke_ul);
+
+            hist.name = "RolkeLL";
+            WriteHistogram(connection, hist, rolke_ll);
+
+            if (*ib=="Energy")
+            {
+                hist.name = "IntegralRolkeUL";
+                WriteHistogram(connection, hist, rolke_int);
+            }
 #endif
         }