pred_err.C
Go to the documentation of this file.
6 #include "CAFAna/Core/Spectrum.h"
7 
9 
10 #include "TFile.h"
11 #include "TH1.h"
12 #include "TCanvas.h"
13 #include "TGaxis.h"
14 #include "TLatex.h"
15 #include "TLegend.h"
16 #include "TLine.h"
17 #include "TSystem.h"
18 #include "TGraphAsymmErrors.h"
19 #include "TLegendEntry.h"
20 #include "TPaveText.h"
21 #include <iostream>
22 #include <iomanip>
23 
24 using namespace ana;
25 
26 
28 {
29  const std::string extrap = beam == "fhc" ? "combo" : "prop";
30  const double pot = beam == "fhc" ? kAna2020FHCPOT : kAna2020RHCPOT;
31  const double livetime = beam == "fhc" ? kAna2020FHCLivetime : kAna2020RHCLivetime;
32 
33  // 2020 Best Fit
34  std::string calcfile = "/nova/ana/3flavor/Ana2020/Fits/BestFits/bestfits_joint_realData_both_systs.root";
35  TFile *calcf = new TFile(calcfile.c_str());
36  auto calcosc = LoadFrom<osc::IOscCalcAdjustable>(calcf, "NH_UO_calc").release();
37 
38  // Load nue2020 prediction
39  auto pred = GetNuePrediction2020(extrap,calcosc,true,beam,false);
40 
41  // Load cosmic prediction
42  TH1* hData_out = GetNueCosmics2020(beam).first->ToTH1(livetime, kLivetime);
43  double cosmic_int = hData_out->Integral();
44 
45  // Load 2020 Systs
47 
48  // 2020 syst shift fit
49  SystShifts * systshifts = LoadFromFile<SystShifts>(calcfile, "NH_UO_systs").release();
50 
51  // Need nue and beam specific systs only.
52  // Also leave empty if we are plotting the error band with the nominal
53  SystShifts myshifts;
54  for(const ISyst* syst: systs)
55  myshifts.SetShift(syst, systshifts->GetShift(syst));
56 
57  // Predict components
58  std::pair<double,double> hSig = {pred->PredictComponentSyst(calcosc, myshifts,ana::Flavors::kNuMuToNuE,ana::Current::kCC,
59  beam=="fhc" ? ana::Sign::kNu : ana::Sign::kAntiNu).ToTH1(pot)->Integral(),
60  pred->PredictComponent(calcosc, ana::Flavors::kNuMuToNuE,ana::Current::kCC,
61  beam=="fhc" ? ana::Sign::kNu : ana::Sign::kAntiNu).ToTH1(pot)->Integral()};
62  std::pair<double,double> hTot = {pred->PredictSyst(calcosc, myshifts).ToTH1(pot)->Integral()+cosmic_int,
63  pred->Predict(calcosc).ToTH1(pot)->Integral()+cosmic_int};
64  std::pair<double,double> hBkg = {hTot.first-hSig.first, hTot.second-hSig.second};
65 
66  // Syst universes
67  std::vector<double> sigups, sigdns;
68  std::vector<double> totups, totdns;
69  std::vector<double> bkgups, bkgdns;
70  for(const ISyst* syst: systs){
71  SystShifts shifts;
72  // +1
73  shifts.SetShift(syst, +1);
74  double htotup = pred->PredictSyst(calcosc, shifts).ToTH1(pot)->Integral() + cosmic_int;
75  double hsigup = pred->PredictComponentSyst(calcosc, shifts, ana::Flavors::kNuMuToNuE, ana::Current::kCC,
76  beam=="fhc" ? ana::Sign::kNu : ana::Sign::kAntiNu).ToTH1(pot)->Integral();
77  double hbkgup = htotup-hsigup;
78  sigups.push_back(hsigup);
79  totups.push_back(htotup);
80  bkgups.push_back(hbkgup);
81 
82  // -1
83  shifts.SetShift(syst, -1);
84  double htotdn = pred->PredictSyst(calcosc, shifts).ToTH1(pot)->Integral()+cosmic_int;
85  double hsigdn = pred->PredictComponentSyst(calcosc, shifts, ana::Flavors::kNuMuToNuE, ana::Current::kCC,
86  beam=="fhc" ? ana::Sign::kNu : ana::Sign::kAntiNu).ToTH1(pot)->Integral();
87  double hbkgdn = htotdn-hsigdn;
88  sigdns.push_back(hsigdn);
89  totdns.push_back(htotdn);
90  bkgdns.push_back(hbkgdn);
91  }
92 
93  totups.push_back(hTot.second+1.143/sqrt(40));
94  totdns.push_back(hTot.second-1.143/sqrt(40));
95  bkgups.push_back(hBkg.second+1.143/sqrt(40));
96  bkgdns.push_back(hBkg.second-1.143/sqrt(40));
97 
98  // Signal
99  double errSigUp = 0;
100  double errSigDn = 0;
101  for(unsigned int systIdx = 0; systIdx < sigups.size(); ++systIdx){
102  double hi = sigups[systIdx]-hSig.second;
103  double lo = sigdns[systIdx]-hSig.second;
104 
105  double min = std::min(hi,lo);
106  double max = std::max(hi,lo);
107  if(max < 0) max=0;
108  if(min > 0) min=0;
109  errSigUp += max*max;
110  errSigDn += min*min;
111  } // end for systIdx
112 
113  // Total
114  double errTotUp = 0;
115  double errTotDn = 0;
116  for(unsigned int systIdx = 0; systIdx < totups.size(); ++systIdx){
117  double hi = totups[systIdx]-hTot.second;
118  double lo = totdns[systIdx]-hTot.second;
119 
120  double min = std::min(hi,lo);
121  double max = std::max(hi,lo);
122  if(max < 0) max=0;
123  if(min > 0) min=0;
124  errTotUp += max*max;
125  errTotDn += min*min;
126  } // end for systIdx
127 
128  // Background
129  double errBkgUp = 0;
130  double errBkgDn = 0;
131  for(unsigned int systIdx = 0; systIdx < bkgups.size(); ++systIdx){
132  double hi = bkgups[systIdx]-hBkg.second;
133  double lo = bkgdns[systIdx]-hBkg.second;
134 
135  double min = std::min(hi,lo);
136  double max = std::max(hi,lo);
137  if(max < 0) max=0;
138  if(min > 0) min=0;
139  errBkgUp += max*max;
140  errBkgDn += min*min;
141  } // end for systIdx
142 
143 
144  std::cout<<hSig.first<<"\t"<<sqrt(errSigUp)<<"\t"<<sqrt(errSigDn)<<std::endl;
145  std::cout<<hTot.first<<"\t"<<sqrt(errTotUp)<<"\t"<<sqrt(errTotDn)<<std::endl;
146  std::cout<<hBkg.first<<"\t"<<sqrt(errBkgUp)<<"\t"<<sqrt(errBkgDn)<<std::endl;
147 
148 
149  // with syst:
150  // 9.1084
151 } // end
T max(const caf::Proxy< T > &a, T b)
TSpline3 lo("lo", xlo, ylo, 12,"0")
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
Antineutrinos-only.
Definition: IPrediction.h:50
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
(&#39; appearance&#39;)
Definition: IPrediction.h:18
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
T sqrt(T number)
Definition: d0nt_math.hpp:156
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
std::pair< Spectrum *, double > GetNueCosmics2020(std::string beam, bool GetFromUPS=false, bool NERSC=false)
Charged-current interactions.
Definition: IPrediction.h:39
TSpline3 hi("hi", xhi, yhi, 18,"0")
#define pot
const IPrediction * GetNuePrediction2020(std::string decomp, osc::IOscCalc *calc, bool corrSysts, std::string beam, bool isFakeData, bool GetFromUPS=false, bool minimizeMemory=true, bool NERSC=false)
const double kAna2020FHCLivetime
Definition: Exposures.h:237
T GetShift(const ISyst *syst) const
const double kAna2020FHCPOT
Definition: Exposures.h:233
void pred_err(const std::string beam)
Definition: pred_err.C:27
const double kAna2020RHCPOT
Definition: Exposures.h:235
Neutrinos-only.
Definition: IPrediction.h:49
OStream cout
Definition: OStream.cxx:6
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
double livetime
Definition: saveFDMCHists.C:21
const double kAna2020RHCLivetime
Definition: Exposures.h:238
std::vector< const ISyst * > get3FlavorAna2020AllSysts(const EAnaType2020 ana, const bool smallgenie, const BeamType2020 beam, const bool isFit, const bool ptExtrap)
T min(const caf::Proxy< T > &a, T b)
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:80
enum BeamMode string