Index: /trunk/FACT++/src/spectrum.cc =================================================================== --- /trunk/FACT++/src/spectrum.cc (revision 19897) +++ /trunk/FACT++/src/spectrum.cc (revision 19898) @@ -1376,5 +1376,5 @@ 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,5 +1822,5 @@ << 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,8 +1831,8 @@ 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,4 +1888,5 @@ map rolke_ul; map rolke_ll; + map rolke_int; for (auto ir=res13.cbegin(); ir!=res13.cend(); ir++) @@ -1894,23 +1895,32 @@ 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,69 +1934,104 @@ } - 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 = "/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 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 = "/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 }