Changeset 19898


Ignore:
Timestamp:
Dec 14, 2019, 11:42:10 PM (4 months ago)
Author:
tbretz
Message:
Updated to include Integral and Integrated Results.
File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/FACT++/src/spectrum.cc

    r19895 r19898  
    13761376        query9.template_defaults["binning"]  = *ib=="theta" ? "BinningTheta" : ("BinningEnergy_"+*ib).c_str();
    13771377        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))";
     1378        query9.template_defaults["binwidth"] = *ib=="theta" ? "1"            : "(POW(10,hi)-POW(10,lo))/1000";
    13791379
    13801380        if (print_queries)
     
    18221822            << indent(3) << spectrum_sql << indent(0) <<
    18231823            ")";
    1824 */
     1824
    18251825        // [ Sa^2/a^2 + Sb^2/b^2 ] * a^2/b^2
    18261826        // [ (sc^2)/c^2+(sb^2)/b^2+(sa^2)/a^2 ]  * a^2*b^2/c^2
     
    18311831
    18321832        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";
     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";
    18371837        if (*ib=="Theta")
    18381838        {
     
    18881888            map<size_t, double> rolke_ul;
    18891889            map<size_t, double> rolke_ll;
     1890            map<size_t, double> rolke_int;
    18901891
    18911892            for (auto ir=res13.cbegin(); ir!=res13.cend(); ir++)
     
    18941895                const mysqlpp::Row &row = *ir;
    18951896
    1896                 const size_t bin     = row[*ib=="Theta" ? ".theta" : ".sparse_est"];
    1897                 const double flux    = row["Flux"];
    1898                 const double error   = row["ErrFlux"];
     1897                const size_t bin      = row[*ib=="Theta" ? ".theta" : ".sparse_est"];
     1898                const double flux     = row["Flux"];
     1899                const double error    = row["ErrFlux"];
    18991900
    19001901#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"];
     1902                const double dat_sig  = row["Signal"];
     1903                const double dat_bg   = row["Background"];
     1904
     1905                const double dat_isig = row["IntegralSignal"];
     1906                const double dat_ibg  = row["IntegralBackground"];
     1907
     1908                const double eff      = row["Efficiency"];
     1909                const double ieff     = row["IntegralEfficiency"];
     1910
     1911                const double areatime = row["AreaTime"];
     1912                const double width    = row["Width"];
    19061913
    19071914                fc.SetMuMax(10*(dat_sig+dat_bg/5)); //has to be higher than sig+bg!!
    1908                 rolke.SetPoissonBkgKnownEff(dat_sig, dat_bg, 5, eff);
    19091915
    19101916                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;
     1917                    feldman_ul[bin] = fc.CalculateUpperLimit(dat_sig, dat_bg/5)/width/areatime/eff;
     1918
     1919                rolke.SetPoissonBkgKnownEff(dat_sig,  dat_bg, 5, eff);
     1920                rolke_ll[bin]  = rolke.GetLowerLimit()/width/areatime;
     1921                rolke_ul[bin]  = rolke.GetUpperLimit()/width/areatime;
     1922
     1923                rolke.SetPoissonBkgKnownEff(dat_isig, dat_ibg, 5, ieff);
     1924                rolke_int[bin] = rolke.GetUpperLimit()/areatime;
    19151925#endif
    19161926                if (verbose>0)
     
    19241934            }
    19251935
    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);
     1936            Histogram hist;
     1937            hist.dir      = "Data/"+*ib;
     1938            hist.table    = table;
     1939            hist.binningx = *ib=="Theta" ? binning_theta : binning_sparse;
     1940            hist.x        = *ib=="Theta" ? ".theta" : ".sparse_est";
     1941            hist.axisx    = *ib=="Theta" ? "Zenith Distance #theta [#circ]" : "Energy lg(E/GeV)";
     1942            hist.stats    = false;
     1943
     1944            const vector<string> types { "", "Integral" };
     1945            for (auto it=types.cbegin(); it!=types.cend(); it++)
     1946            {
     1947                hist.axisy = "Counts";
     1948                if (*ib=="Energy")
     1949                    hist.axisy += " (E>E_{lo})";
     1950
     1951                hist.title = "";
     1952                hist.v     = *it+"Signal";
     1953                hist.err   = "Err"+*it+"Signal";
     1954                WriteHistogram(connection, hist);
     1955
     1956                hist.title = "";
     1957                hist.v     = *it+"Background";
     1958                hist.err   = "Err"+*it+"Background";
     1959                WriteHistogram(connection, hist);
     1960
     1961                hist.title = "";
     1962                hist.v     = *it+"Excess";
     1963                hist.err   = "Err"+*it+"Excess";
     1964                WriteHistogram(connection, hist);
     1965
     1966                hist.title = "";
     1967                hist.v     = *it+"Significance";
     1968                hist.err   = "";
     1969                hist.axisy = "#sigma";
     1970                if (*ib=="Energy")
     1971                    hist.axisy += " (E>E_{lo})";
     1972                WriteHistogram(connection, hist);
     1973
     1974                hist.title = "";
     1975                hist.v     = *it+"AvgEnergyEst";
     1976                hist.err   = "";
     1977                hist.axisy = "<E_{est}>/GeV";
     1978                if (*ib=="Energy")
     1979                    hist.axisy += " (E>E_{lo})";
     1980                WriteHistogram(connection, hist);
     1981
     1982                hist.title = "";
     1983                hist.v     = *it+"ExcessRatio";
     1984                hist.err   = "Err"+*it+"ExcessRatio";
     1985                hist.axisy = "Ratio";
     1986                if (*ib=="Energy")
     1987                    hist.axisy += " (E>E_{lo})";
     1988                WriteHistogram(connection, hist);
     1989            }
     1990
     1991            hist.axisy = "dN/dE ";
     1992            if (*ib=="Energy")
     1993                hist.axisy += "(E>E_{lo}) ";
     1994            hist.axisy += *ib=="Theta" ? "[cm^{-2} s^{-1}]" : "[cm^{-2} s^{-1} TeV^{-1}]";
     1995
     1996            hist.name = "Spectrum";
     1997            hist.v    = "Flux";
     1998            hist.err  = "ErrFlux";
     1999            WriteHistogram(connection, hist);
     2000
     2001            if (*ib=="Energy")
     2002            {
     2003                hist.axisy = "dN/dE (E>E_{lo}) [cm^{-2} s^{-1}]";
     2004
     2005                hist.name  = "IntegralSpectrum";
     2006                hist.v     = "IntegralFlux";
     2007                hist.err   = "ErrIntegralFlux";
     2008                WriteHistogram(connection, hist);
     2009
     2010                hist.name  = "IntegratedSpectrum";
     2011                hist.v     = "IntegratedFlux";
     2012                hist.err   = "ErrIntegratedFlux";
     2013                WriteHistogram(connection, hist);
     2014            }
    19762015
    19772016#ifdef HAVE_ROOT
    1978             hist_zd.axisy    = *ib=="Theta" ? "UL [cm^{-2} s^{-1}]" : "UL [cm^{-2} s^{-1} TeV^{-1}]";
     2017            hist.axisy = *ib=="Theta" ? "UL [cm^{-2} s^{-1}]" : "UL [cm^{-2} s^{-1} TeV^{-1}]";
    19792018
    19802019            if (feldman)
    19812020            {
    1982                 hist_zd.name = "FeldmanCousins";
    1983                 WriteHistogram(connection, hist_zd, feldman_ul);
     2021                hist.name = "FeldmanCousins";
     2022                WriteHistogram(connection, hist, feldman_ul);
    19842023            }
    19852024
    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);
     2025            hist.name = "RolkeUL";
     2026            WriteHistogram(connection, hist, rolke_ul);
     2027
     2028            hist.name = "RolkeLL";
     2029            WriteHistogram(connection, hist, rolke_ll);
     2030
     2031            if (*ib=="Energy")
     2032            {
     2033                hist.name = "IntegralRolkeUL";
     2034                WriteHistogram(connection, hist, rolke_int);
     2035            }
    19912036#endif
    19922037        }
Note: See TracChangeset for help on using the changeset viewer.