FluxMultiverseSyst.cxx
Go to the documentation of this file.
2 
3 #include <iostream>
4 #include <string>
5 
6 #include "CAFAna/Core/Spectrum.h"
7 #include "CAFAna/Systs/Systs.h"
8 
10 
11 #include "TH1.h"
12 #include "TObjString.h"
13 
14 using namespace std;
15 
16 namespace ana
17 {
18  //////
19  FluxMultiverseSyst::FluxMultiverseSyst(SpectrumLoaderBase& loader,
20  const HistAxis& axis,
21  const Cut& cut,
22  const SystShifts& shift,
23  const Var& cv_wei,
24  const std::vector<Var>& vuniv_wei):
25  fNUniverses(0),fUpperSigma(NULL),fLowerSigma(NULL),fBoundariesFound(false),
26  fNormUpperSigma(NULL), fNormLowerSigma(NULL){
27 
28  for(unsigned int i = 0; i < vuniv_wei.size(); i++){
29  fSpectra.push_back( new Spectrum(loader, axis, cut, shift, vuniv_wei[i]));
30  }
31  fNUniverses = vuniv_wei.size();
32  // last position is the cv
33  fSpectra.push_back( new Spectrum(loader, axis, cut, shift, cv_wei) );
34  }
35 
36  //////
38  std::vector<std::unique_ptr<ana::Spectrum>>& spectra, bool fromfile): fLoadFromFile(fromfile)
39  {
41  for(auto& it: spectra){
42  fSpectra.push_back(it.release());
43  }
44 
45  }
46 
47  FluxMultiverseSyst::FluxMultiverseSyst(std::vector<std::unique_ptr<ana::Spectrum>>& spectra){
48  fNUniverses = spectra.size()-1;
49  for(auto& it: spectra){
50  fSpectra.push_back(it.release());
51  }
52  }
53 
55  if(!fLoadFromFile){
56  for (auto it = fSpectra.begin(); it != fSpectra.end(); ++it) delete *it;
57  }
58  delete fUpperSigma;
59  delete fLowerSigma;
60  delete fNormUpperSigma;
61  delete fNormLowerSigma;
62  }
63 
65  return fSpectra[fNUniverses];
66  }
67 
69  {
71  return fNormUpperSigma;
72  }
73 
75  {
77  return fNormLowerSigma;
78  }
79 
80  TGraphAsymmErrors* FluxMultiverseSyst::ToAreaNormalizedTH1(double pot, int col, int errCol,
81  float headroom, bool newaxis){
82  if(col == -1){
83  col = kRed;
84  errCol = kRed-10;
85  }
86  else if(errCol == -1) errCol = col-7; // hopefully a lighter version
87 
88  TH1* cv = CV()->ToTH1(pot);
89 
90  std::vector<TH1*> ups, dns;
91  ups.push_back(fUpperSigma->ToTH1(pot));
92  dns.push_back(fLowerSigma->ToTH1(pot));
93 
94  cv->SetLineColor(col);
95  cv->GetXaxis()->CenterTitle();
96  cv->GetYaxis()->CenterTitle();
97 
98  double yMax = cv->GetBinContent(cv->GetMaximumBin());
99 
100  TGraphAsymmErrors* g = new TGraphAsymmErrors;
101 
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);
105 
106  const double w = cv->GetXaxis()->GetBinWidth(binIdx);
107 
108  double errUp = 0;
109  double errDn = 0;
110 
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);
114 
115  if(lo <= 0 && hi <= 0) std::swap(lo, hi);
116 
117  errUp += hi*hi;
118  errDn += lo*lo;
119 
120  // TODO: what happens if they're both high or both low?
121  } // end for systIdx
122 
123 
124  g->SetPointError(binIdx, w/2, w/2, sqrt(errDn), sqrt(errUp));
125  } // end for i
126 
127  g->SetFillColor(errCol);
128  //~ g->Draw("e2 same");
129  g->GetYaxis()->SetRangeUser(0, headroom*yMax);
130  cv->GetYaxis()->SetRangeUser(0, headroom*yMax);
131 
132  //~ nom->Draw("hist ][ same");
133  return g;
134  }
135 
137  {
138  return ToTH1(CV()->POT());
139  }
140 
141  // return a TH1* with Sigma of all universes as the error band
143  TH1* hCV = CV()->ToTH1(pot);
144  vector<TH1*> hUniverses;
145  map< int,float > cv_counts_in_bins;
146  map< int, vector<float> > univ_counts_in_bins;
147 
148  unsigned int count = 0;
149  for(auto& it: fSpectra){
150 
151  if(count!=fNUniverses){
152  TH1* curhist = it->ToTH1(pot);
153  hUniverses.push_back(curhist);
154  int nbinsx = curhist->GetNbinsX();
155  for(int i = 1; i <= nbinsx; i++){
156  univ_counts_in_bins[i].push_back(curhist->GetBinContent(i));
157  }
158  }
159  else if(count==fNUniverses){
160  int nbinsx = hCV->GetNbinsX();
161  for(int i = 1; i <= nbinsx; i++){
162  cv_counts_in_bins[i] = hCV->GetBinContent(i);
163  }
164  }
165  count++;
166  } //end loop on spectra.
167 
168  for(auto& it: univ_counts_in_bins){
169  float cv_val = hCV->GetBinContent( it.first);
170  if(cv_val <=0)continue;
171  float err = 0.0;
172  for(unsigned int iuniv=0;iuniv<fNUniverses;iuniv++){
173  err += pow(it.second[iuniv] - cv_val,2);
174  }
175  err /= float(fNUniverses);
176  err = sqrt(err);
177  hCV->SetBinError(it.first, err);
178  }
179  hCV->SetLineColor(kRed);
180  hCV->SetMarkerColor(kRed);
181 
182  return hCV;
183  }
184 
185  // return upper boundary
187  {
189  return fUpperSigma;
190  }
191 
192  // function to fill the boudary spactra encircling the error band
194  Eigen::ArrayXd aCV = CV()->GetEigen(CV()->POT());
195  map< int,float > cv_events_in_bins;
196  map<int, vector<float> > univ_events_in_bins;
197  unsigned int count = 0;
198  for(auto& it: fSpectra){
199  if(count!=fNUniverses){
200  Eigen::ArrayXd curhist = it->GetEigen(it->POT());
201  int nbinsx = curhist.size()-2;
202  for(int i = 1; i <= nbinsx; i++){
203  univ_events_in_bins[i].push_back(curhist[i]);
204  }
205  }
206  else if(count==fNUniverses){
207  int nbinsx = aCV.size()-2;
208  for(int i = 1; i <= nbinsx; i++){
209  cv_events_in_bins[i] = aCV[i];
210  }
211  }
212  count++;
213  }
214 
215  Eigen::ArrayXd aUp = aCV;
216  Eigen::ArrayXd aLow = aCV;
217  int binIt = 1;
218  for(auto& it: univ_events_in_bins){
219  double cv_val = aCV[it.first];
220  double err = 0.0;
221  for(unsigned int iuniv=0;iuniv<fNUniverses;iuniv++){
222  err += pow(it.second[iuniv] - cv_val,2);
223  }
224  err /= double(fNUniverses);
225  err = sqrt(err);
226  aUp[binIt] = cv_val+err;
227  aLow[binIt] = cv_val-err;
228  binIt++;
229  }
230 
231  // release memory if they are already allocated
232  delete fUpperSigma;
233  delete fLowerSigma;
234  delete fNormUpperSigma;
235  delete fNormLowerSigma;
236  fUpperSigma = new Spectrum(std::move(aUp),
237  HistAxis(CV()->GetLabels(), CV()->GetBinnings()),
238  fSpectra[0]->POT(), fSpectra[0]->Livetime());
239  fLowerSigma = new Spectrum(std::move(aLow),
240  HistAxis(CV()->GetLabels(), CV()->GetBinnings()),
241  fSpectra[0]->POT(), fSpectra[0]->Livetime());
242  fBoundariesFound = true;
243  Spectrum cv = *CV();
244  Spectrum tempMC = cv;
249  }
250 
251  // return lower boundary
253  {
255  return fLowerSigma;
256  }
257 
258  //----------------------------------------------------------------------
259  void FluxMultiverseSyst::SaveTo(TDirectory* dir, const std::string& name) const {
260  TDirectory *tmp = gDirectory;
261 
262  dir = dir->mkdir(name.c_str()); // switch to subdir
263  dir->cd();
264 
265  TObjString universes(Form("%lu",fSpectra.size()));
266  universes.Write("nUniverses");
267 
268  unsigned int count = 0;
269  for(auto& it: fSpectra){
270  if(count != fNUniverses){
271  it->SaveTo(dir, "flux_universe_"+std::to_string(count));
272  }
273  if(count == fNUniverses){
274  it->SaveTo(dir, "flux_cv");
275  }
276  count++;
277  }// end loop over spectra
278 
279  dir->Write();
280  delete dir;
281 
282  tmp->cd();
283  }// end of SaveTo
284 
285 
286  std::unique_ptr<FluxMultiverseSyst> FluxMultiverseSyst::LoadFrom(TDirectory* dir, const std::string& name) {
287  dir = dir->GetDirectory(name.c_str()); // switch to subdir
288  assert(dir);
289 
290  TObjString* universes = (TObjString*)dir->Get("nUniverses");
291  unsigned int nUniverses = universes->GetString().Atoi();
292 
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)));
296  }
297  spectra.push_back(ana::LoadFrom<Spectrum>(dir, "flux_cv"));
298 
299  delete dir;
300 
301  return std::make_unique<FluxMultiverseSyst>(nUniverses, spectra, true);
302 
303  }// end of LoadFrom
304 }
Spectrum * fNormUpperSigma
spectra of area normalized upper rms and lower rms
const XML_Char * name
Definition: expat.h:151
static std::unique_ptr< FluxMultiverseSyst > LoadFrom(TDirectory *dir, const std::string &name)
TSpline3 lo("lo", xlo, ylo, 12,"0")
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
set< int >::iterator it
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.
Definition: Spectrum.cxx:148
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
void OverridePOT(double newpot)
DO NOT USE UNLESS YOU ARE 110% CERTAIN THERE ISN&#39;T A BETTER WAY!
Definition: Spectrum.h:225
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
void SaveTo(TDirectory *dir, const std::string &name) const
const Eigen::ArrayXd & GetEigen() const
NB these don&#39;t have POT scaling. For expert high performance ops only!
Definition: Spectrum.h:188
Float_t tmp
Definition: plot.C:36
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:272
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.
Definition: Spectrum.h:40
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
unsigned int fNUniverses
number of universes to generate
TSpline3 hi("hi", xhi, yhi, 18,"0")
Int_t col[ntarg]
Definition: Style.C:29
const std::string cv[Ncv]
#define pot
TGraphAsymmErrors * ToAreaNormalizedTH1(double pot, int col=-1, int errCol=-1, float headroom=1.3, bool newaxis=true)
loader
Definition: demo0.py:10
std::vector< float > Spectrum
Definition: Constants.h:570
double POT() const
Definition: Spectrum.h:219
std::vector< double > POT
Base class for the various types of spectrum loader.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const Cut cut
Definition: exporter_fd.C:30
Int_t nbinsx
Definition: plot.C:23
TDirectory * dir
Definition: macro.C:5
assert(nhit_max >=nhit_nbins)
Spectrum * NormUpperSigma()
return area normalized sigma extremes
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
Spectrum * fUpperSigma
spectra of upper sigma and lower sigma
const int nuniverses
Spectrum * UpperSigma()
return sigma extremes
Float_t w
Definition: plot.C:20