Public Member Functions | Static Public Member Functions | Private Member Functions | Private Attributes | List of all members
ana::FluxMultiverseSyst Class Reference

Obtain PPFX GENIE systematic band. More...

#include "/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N20-11-30/CAFAna/XSec/FluxMultiverseSyst.h"

Public Member Functions

 FluxMultiverseSyst ()
 
 FluxMultiverseSyst (SpectrumLoaderBase &loader, const HistAxis &axis, const Cut &cut, const SystShifts &shift, const Var &cv_wei=kUnweighted, const std::vector< Var > &vuniv_wei=std::vector< Var >())
 
 FluxMultiverseSyst (unsigned int nuniverses, std::vector< std::unique_ptr< ana::Spectrum >> &spectra, bool fromfile=false)
 
 FluxMultiverseSyst (std::vector< std::unique_ptr< ana::Spectrum >> &spectra)
 
 ~FluxMultiverseSyst ()
 
TH1 * ToTH1 ()
 
TH1 * ToTH1 (double pot)
 
TGraphAsymmErrors * ToAreaNormalizedTH1 (double pot, int col=-1, int errCol=-1, float headroom=1.3, bool newaxis=true)
 
std::vector< Spectrum * > AllUniverses ()
 return all universes More...
 
SpectrumUpperSigma ()
 return sigma extremes More...
 
SpectrumLowerSigma ()
 
SpectrumCV ()
 
SpectrumNormUpperSigma ()
 return area normalized sigma extremes More...
 
SpectrumNormLowerSigma ()
 
void SaveTo (TDirectory *dir, const std::string &name) const
 

Static Public Member Functions

static std::unique_ptr< FluxMultiverseSystLoadFrom (TDirectory *dir, const std::string &name)
 

Private Member Functions

void FindSigmaBoundaries ()
 fUpperSigma and fLowerSigma More...
 

Private Attributes

unsigned int fNUniverses
 number of universes to generate More...
 
bool fLoadFromFile = false
 
std::vector< Spectrum * > fSpectra
 spectra of all universes including the central value: More...
 
SpectrumfUpperSigma = NULL
 spectra of upper sigma and lower sigma More...
 
SpectrumfLowerSigma = NULL
 
bool fBoundariesFound
 
SpectrumfNormUpperSigma = NULL
 spectra of area normalized upper rms and lower rms More...
 
SpectrumfNormLowerSigma = NULL
 

Detailed Description

Obtain PPFX GENIE systematic band.

Definition at line 19 of file FluxMultiverseSyst.h.

Constructor & Destructor Documentation

ana::FluxMultiverseSyst::FluxMultiverseSyst ( )
inline
ana::FluxMultiverseSyst::FluxMultiverseSyst ( SpectrumLoaderBase loader,
const HistAxis axis,
const Cut cut,
const SystShifts shift,
const Var cv_wei = kUnweighted,
const std::vector< Var > &  vuniv_wei = std::vector<Var>() 
)

Definition at line 19 of file FluxMultiverseSyst.cxx.

References fNUniverses, fSpectra, and MECModelEnuComparisons::i.

24  :
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  }
Spectrum * fNormUpperSigma
spectra of area normalized upper rms and lower rms
std::vector< Spectrum * > fSpectra
spectra of all universes including the central value:
unsigned int fNUniverses
number of universes to generate
loader
Definition: demo0.py:10
std::vector< float > Spectrum
Definition: Constants.h:610
const Cut cut
Definition: exporter_fd.C:30
Spectrum * fUpperSigma
spectra of upper sigma and lower sigma
ana::FluxMultiverseSyst::FluxMultiverseSyst ( unsigned int  nuniverses,
std::vector< std::unique_ptr< ana::Spectrum >> &  spectra,
bool  fromfile = false 
)

Definition at line 37 of file FluxMultiverseSyst.cxx.

References fNUniverses, fSpectra, it, and nuniverses.

39  {
41  for(auto& it: spectra){
42  fSpectra.push_back(it.release());
43  }
44 
45  }
set< int >::iterator it
std::vector< Spectrum * > fSpectra
spectra of all universes including the central value:
unsigned int fNUniverses
number of universes to generate
const int nuniverses
ana::FluxMultiverseSyst::FluxMultiverseSyst ( std::vector< std::unique_ptr< ana::Spectrum >> &  spectra)

Definition at line 47 of file FluxMultiverseSyst.cxx.

References fNUniverses, fSpectra, and it.

47  {
48  fNUniverses = spectra.size()-1;
49  for(auto& it: spectra){
50  fSpectra.push_back(it.release());
51  }
52  }
set< int >::iterator it
std::vector< Spectrum * > fSpectra
spectra of all universes including the central value:
unsigned int fNUniverses
number of universes to generate
ana::FluxMultiverseSyst::~FluxMultiverseSyst ( )

Definition at line 54 of file FluxMultiverseSyst.cxx.

References fLoadFromFile, fLowerSigma, fNormLowerSigma, fNormUpperSigma, fSpectra, fUpperSigma, and it.

Referenced by FluxMultiverseSyst().

54  {
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  }
Spectrum * fNormUpperSigma
spectra of area normalized upper rms and lower rms
set< int >::iterator it
std::vector< Spectrum * > fSpectra
spectra of all universes including the central value:
Spectrum * fUpperSigma
spectra of upper sigma and lower sigma

Member Function Documentation

std::vector<Spectrum*> ana::FluxMultiverseSyst::AllUniverses ( )
inline

return all universes

Definition at line 45 of file FluxMultiverseSyst.h.

References CV(), dir, FindSigmaBoundaries(), fSpectra, LoadFrom(), LowerSigma(), NormLowerSigma(), NormUpperSigma(), SaveTo(), string, and UpperSigma().

45 {return fSpectra;}
std::vector< Spectrum * > fSpectra
spectra of all universes including the central value:
Spectrum * ana::FluxMultiverseSyst::CV ( )

Definition at line 64 of file FluxMultiverseSyst.cxx.

References fNUniverses, and fSpectra.

Referenced by AllUniverses(), FindSigmaBoundaries(), ToAreaNormalizedTH1(), and ToTH1().

64  {
65  return fSpectra[fNUniverses];
66  }
std::vector< Spectrum * > fSpectra
spectra of all universes including the central value:
unsigned int fNUniverses
number of universes to generate
void ana::FluxMultiverseSyst::FindSigmaBoundaries ( )
private

fUpperSigma and fLowerSigma

Definition at line 193 of file FluxMultiverseSyst.cxx.

References confusionMatrixTree::count, cv, CV(), fBoundariesFound, fLowerSigma, fNormLowerSigma, fNormUpperSigma, fNUniverses, fSpectra, fUpperSigma, ana::Spectrum::GetEigen(), MECModelEnuComparisons::i, ana::Spectrum::Integral(), it, nbinsx, ana::Spectrum::OverridePOT(), POT, ana::Spectrum::POT(), cet::pow(), and std::sqrt().

Referenced by AllUniverses(), LowerSigma(), NormLowerSigma(), NormUpperSigma(), and UpperSigma().

193  {
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;
246  fNormUpperSigma->OverridePOT(fNormUpperSigma->POT() * fNormUpperSigma->Integral(1) / tempMC.Integral(1));
248  fNormLowerSigma->OverridePOT(fNormLowerSigma->POT() * fNormLowerSigma->Integral(1) / tempMC.Integral(1));
249  }
Spectrum * fNormUpperSigma
spectra of area normalized upper rms and lower rms
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
set< int >::iterator it
void OverridePOT(double newpot)
DO NOT USE UNLESS YOU ARE 110% CERTAIN THERE ISN&#39;T A BETTER WAY!
Definition: Spectrum.h:233
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
const Eigen::ArrayXd & GetEigen() const
NB these don&#39;t have POT scaling. For expert high performance ops only!
Definition: Spectrum.h:189
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:249
std::vector< Spectrum * > fSpectra
spectra of all universes including the central value:
unsigned int fNUniverses
number of universes to generate
const std::string cv[Ncv]
std::vector< float > Spectrum
Definition: Constants.h:610
double POT() const
Definition: Spectrum.h:227
std::vector< double > POT
Int_t nbinsx
Definition: plot.C:23
Spectrum * fUpperSigma
spectra of upper sigma and lower sigma
std::unique_ptr< FluxMultiverseSyst > ana::FluxMultiverseSyst::LoadFrom ( TDirectory *  dir,
const std::string name 
)
static

Definition at line 286 of file FluxMultiverseSyst.cxx.

References ana::assert(), dir, MECModelEnuComparisons::i, and art::to_string().

Referenced by AllUniverses().

286  {
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
const XML_Char * name
Definition: expat.h:151
TDirectory * dir
Definition: macro.C:5
assert(nhit_max >=nhit_nbins)
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
Spectrum * ana::FluxMultiverseSyst::LowerSigma ( )

Definition at line 252 of file FluxMultiverseSyst.cxx.

References fBoundariesFound, FindSigmaBoundaries(), and fLowerSigma.

Referenced by AllUniverses().

253  {
255  return fLowerSigma;
256  }
void FindSigmaBoundaries()
fUpperSigma and fLowerSigma
Spectrum * ana::FluxMultiverseSyst::NormLowerSigma ( )

Definition at line 74 of file FluxMultiverseSyst.cxx.

References fBoundariesFound, FindSigmaBoundaries(), and fNormLowerSigma.

Referenced by AllUniverses().

75  {
77  return fNormLowerSigma;
78  }
void FindSigmaBoundaries()
fUpperSigma and fLowerSigma
Spectrum * ana::FluxMultiverseSyst::NormUpperSigma ( )

return area normalized sigma extremes

Definition at line 68 of file FluxMultiverseSyst.cxx.

References fBoundariesFound, FindSigmaBoundaries(), and fNormUpperSigma.

Referenced by AllUniverses().

69  {
71  return fNormUpperSigma;
72  }
Spectrum * fNormUpperSigma
spectra of area normalized upper rms and lower rms
void FindSigmaBoundaries()
fUpperSigma and fLowerSigma
void ana::FluxMultiverseSyst::SaveTo ( TDirectory *  dir,
const std::string name 
) const

Definition at line 259 of file FluxMultiverseSyst.cxx.

References confusionMatrixTree::count, dir, fNUniverses, fSpectra, it, tmp, and art::to_string().

Referenced by AllUniverses().

259  {
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
const XML_Char * name
Definition: expat.h:151
set< int >::iterator it
Float_t tmp
Definition: plot.C:36
std::vector< Spectrum * > fSpectra
spectra of all universes including the central value:
unsigned int fNUniverses
number of universes to generate
TDirectory * dir
Definition: macro.C:5
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
TGraphAsymmErrors * ana::FluxMultiverseSyst::ToAreaNormalizedTH1 ( double  pot,
int  col = -1,
int  errCol = -1,
float  headroom = 1.3,
bool  newaxis = true 
)

Definition at line 80 of file FluxMultiverseSyst.cxx.

References cv, CV(), fLowerSigma, fUpperSigma, MECModelEnuComparisons::g, hi(), kRed, lo(), std::sqrt(), std::swap(), ana::Spectrum::ToTH1(), w, and submit_syst::y.

Referenced by FluxMultiverseSyst().

81  {
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  }
TSpline3 lo("lo", xlo, ylo, 12,"0")
enum BeamMode kRed
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
T sqrt(T number)
Definition: d0nt_math.hpp:156
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
TSpline3 hi("hi", xhi, yhi, 18,"0")
Int_t col[ntarg]
Definition: Style.C:29
const std::string cv[Ncv]
#define pot
Spectrum * fUpperSigma
spectra of upper sigma and lower sigma
Float_t w
Definition: plot.C:20
TH1 * ana::FluxMultiverseSyst::ToTH1 ( )

return the nominal histogram with the RMS of all other universes as error band

Definition at line 136 of file FluxMultiverseSyst.cxx.

References CV(), and POT.

Referenced by FluxMultiverseSyst().

137  {
138  return ToTH1(CV()->POT());
139  }
std::vector< double > POT
TH1 * ana::FluxMultiverseSyst::ToTH1 ( double  pot)

Definition at line 142 of file FluxMultiverseSyst.cxx.

References confusionMatrixTree::count, CV(), check_time_usage::float, fNUniverses, fSpectra, MECModelEnuComparisons::i, it, kRed, nbinsx, cet::pow(), std::sqrt(), and ana::Spectrum::ToTH1().

142  {
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  }
enum BeamMode kRed
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
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
std::vector< Spectrum * > fSpectra
spectra of all universes including the central value:
unsigned int fNUniverses
number of universes to generate
#define pot
Int_t nbinsx
Definition: plot.C:23
Spectrum * ana::FluxMultiverseSyst::UpperSigma ( )

return sigma extremes

Definition at line 186 of file FluxMultiverseSyst.cxx.

References fBoundariesFound, FindSigmaBoundaries(), and fUpperSigma.

Referenced by AllUniverses().

187  {
189  return fUpperSigma;
190  }
void FindSigmaBoundaries()
fUpperSigma and fLowerSigma
Spectrum * fUpperSigma
spectra of upper sigma and lower sigma

Member Data Documentation

bool ana::FluxMultiverseSyst::fBoundariesFound
private
bool ana::FluxMultiverseSyst::fLoadFromFile = false
private

Definition at line 63 of file FluxMultiverseSyst.h.

Referenced by ~FluxMultiverseSyst().

Spectrum* ana::FluxMultiverseSyst::fLowerSigma = NULL
mutableprivate
Spectrum* ana::FluxMultiverseSyst::fNormLowerSigma = NULL
mutableprivate

Definition at line 73 of file FluxMultiverseSyst.h.

Referenced by FindSigmaBoundaries(), NormLowerSigma(), and ~FluxMultiverseSyst().

Spectrum* ana::FluxMultiverseSyst::fNormUpperSigma = NULL
mutableprivate

spectra of area normalized upper rms and lower rms

Definition at line 72 of file FluxMultiverseSyst.h.

Referenced by FindSigmaBoundaries(), NormUpperSigma(), and ~FluxMultiverseSyst().

unsigned int ana::FluxMultiverseSyst::fNUniverses
private

number of universes to generate

Definition at line 61 of file FluxMultiverseSyst.h.

Referenced by CV(), FindSigmaBoundaries(), FluxMultiverseSyst(), SaveTo(), and ToTH1().

std::vector<Spectrum*> ana::FluxMultiverseSyst::fSpectra
private

spectra of all universes including the central value:

Definition at line 66 of file FluxMultiverseSyst.h.

Referenced by AllUniverses(), CV(), FindSigmaBoundaries(), FluxMultiverseSyst(), SaveTo(), ToTH1(), and ~FluxMultiverseSyst().

Spectrum* ana::FluxMultiverseSyst::fUpperSigma = NULL
mutableprivate

spectra of upper sigma and lower sigma

Definition at line 68 of file FluxMultiverseSyst.h.

Referenced by FindSigmaBoundaries(), ToAreaNormalizedTH1(), UpperSigma(), and ~FluxMultiverseSyst().


The documentation for this class was generated from the following files: