12 #include "TObjString.h" 24 const std::vector<Var>& vuniv_wei):
25 fNUniverses(0),fUpperSigma(NULL),fLowerSigma(NULL),fBoundariesFound(false),
26 fNormUpperSigma(NULL), fNormLowerSigma(NULL){
28 for(
unsigned int i = 0;
i < vuniv_wei.size();
i++){
41 for(
auto&
it: spectra){
49 for(
auto&
it: spectra){
81 float headroom,
bool newaxis){
86 else if(errCol == -1) errCol = col-7;
90 std::vector<TH1*> ups, dns;
94 cv->SetLineColor(col);
95 cv->GetXaxis()->CenterTitle();
96 cv->GetYaxis()->CenterTitle();
98 double yMax = cv->GetBinContent(cv->GetMaximumBin());
100 TGraphAsymmErrors*
g =
new TGraphAsymmErrors;
102 for(
int binIdx = 0; binIdx < cv->GetNbinsX()+2; ++binIdx){
103 const double y = cv->GetBinContent(binIdx);
104 g->SetPoint(binIdx, cv->GetXaxis()->GetBinCenter(binIdx),
y);
106 const double w = cv->GetXaxis()->GetBinWidth(binIdx);
111 for(
unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
112 double hi = ups[systIdx]->GetBinContent(binIdx)-
y;
113 double lo = y-dns[systIdx]->GetBinContent(binIdx);
115 if(lo <= 0 && hi <= 0)
std::swap(lo, hi);
124 g->SetPointError(binIdx, w/2, w/2,
sqrt(errDn),
sqrt(errUp));
127 g->SetFillColor(errCol);
129 g->GetYaxis()->SetRangeUser(0, headroom*yMax);
130 cv->GetYaxis()->SetRangeUser(0, headroom*yMax);
144 vector<TH1*> hUniverses;
145 map< int,float > cv_counts_in_bins;
146 map< int, vector<float> > univ_counts_in_bins;
148 unsigned int count = 0;
152 TH1* curhist =
it->ToTH1(pot);
153 hUniverses.push_back(curhist);
154 int nbinsx = curhist->GetNbinsX();
156 univ_counts_in_bins[
i].push_back(curhist->GetBinContent(
i));
160 int nbinsx = hCV->GetNbinsX();
162 cv_counts_in_bins[
i] = hCV->GetBinContent(
i);
168 for(
auto&
it: univ_counts_in_bins){
169 float cv_val = hCV->GetBinContent(
it.first);
170 if(cv_val <=0)
continue;
172 for(
unsigned int iuniv=0;iuniv<
fNUniverses;iuniv++){
173 err +=
pow(
it.second[iuniv] - cv_val,2);
175 err /=
float(fNUniverses);
177 hCV->SetBinError(
it.first, err);
179 hCV->SetLineColor(
kRed);
180 hCV->SetMarkerColor(
kRed);
195 map< int,float > cv_events_in_bins;
196 map<int, vector<float> > univ_events_in_bins;
197 unsigned int count = 0;
200 Eigen::ArrayXd curhist =
it->GetEigen(
it->POT());
201 int nbinsx = curhist.size()-2;
203 univ_events_in_bins[
i].push_back(curhist[
i]);
207 int nbinsx = aCV.size()-2;
209 cv_events_in_bins[
i] = aCV[
i];
215 Eigen::ArrayXd aUp = aCV;
216 Eigen::ArrayXd aLow = aCV;
218 for(
auto&
it: univ_events_in_bins){
219 double cv_val = aCV[
it.first];
221 for(
unsigned int iuniv=0;iuniv<
fNUniverses;iuniv++){
222 err +=
pow(
it.second[iuniv] - cv_val,2);
224 err /= double(fNUniverses);
226 aUp[binIt] = cv_val+err;
227 aLow[binIt] = cv_val-err;
238 fSpectra[0]->
POT(), fSpectra[0]->Livetime());
241 fSpectra[0]->
POT(), fSpectra[0]->Livetime());
260 TDirectory *
tmp = gDirectory;
262 dir = dir->mkdir(name.c_str());
265 TObjString universes(Form(
"%lu",
fSpectra.size()));
266 universes.Write(
"nUniverses");
268 unsigned int count = 0;
274 it->SaveTo(dir,
"flux_cv");
287 dir = dir->GetDirectory(name.c_str());
290 TObjString* universes = (TObjString*)dir->Get(
"nUniverses");
291 unsigned int nUniverses = universes->GetString().Atoi();
293 std::vector<std::unique_ptr<ana::Spectrum>> spectra;
294 for(
unsigned int i = 0;
i < nUniverses;
i++){
295 spectra.push_back(ana::LoadFrom<Spectrum>(dir,
"flux_universe_"+
std::to_string(
i)));
297 spectra.push_back(ana::LoadFrom<Spectrum>(dir,
"flux_cv"));
301 return std::make_unique<FluxMultiverseSyst>(nUniverses, spectra,
true);
Spectrum * fNormUpperSigma
spectra of area normalized upper rms and lower rms
static std::unique_ptr< FluxMultiverseSyst > LoadFrom(TDirectory *dir, const std::string &name)
TSpline3 lo("lo", xlo, ylo, 12,"0")
_HistAxis< Var > HistAxis
Cuts and Vars for the 2020 FD DiF Study.
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Simple record of shifts applied to systematic parameters.
void OverridePOT(double newpot)
DO NOT USE UNLESS YOU ARE 110% CERTAIN THERE ISN'T A BETTER WAY!
void SaveTo(TDirectory *dir, const std::string &name) const
const Eigen::ArrayXd & GetEigen() const
NB these don't have POT scaling. For expert high performance ops only!
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
void FindSigmaBoundaries()
fUpperSigma and fLowerSigma
std::vector< Spectrum * > fSpectra
spectra of all universes including the central value:
Representation of a spectrum in any variable, with associated POT.
Spectrum * NormLowerSigma()
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
unsigned int fNUniverses
number of universes to generate
TSpline3 hi("hi", xhi, yhi, 18,"0")
Spectrum * fNormLowerSigma
const std::string cv[Ncv]
TGraphAsymmErrors * ToAreaNormalizedTH1(double pot, int col=-1, int errCol=-1, float headroom=1.3, bool newaxis=true)
std::vector< float > Spectrum
std::vector< double > POT
Base class for the various types of spectrum loader.
assert(nhit_max >=nhit_nbins)
Spectrum * NormUpperSigma()
return area normalized sigma extremes
std::string to_string(ModuleType mt)
Spectrum * fUpperSigma
spectra of upper sigma and lower sigma
Spectrum * UpperSigma()
return sigma extremes