CalcRWithSysts.C
Go to the documentation of this file.
1 // Macro to plot FD Data/MC spectra with systematic error band
2 
5 #include "CAFAna/Core/Spectrum.h"
8 #include "CAFAna/Systs/Systs.h"
9 
10 #include "OscLib/OscCalcSterile.h"
11 
12 #include "TFile.h"
13 #include "TH1.h"
14 
15 #include <iostream>
16 #include <string>
17 #include <vector>
18 
19 using namespace ana;
20 
21 double CalcR(const Spectrum& sData, const Spectrum& sMCCC, const Spectrum& sMCNC);
22 double CalcDiffR(double r, double rdiff);
23 
25 {
26  TH1::AddDirectory(0);
27 
28  // Create vector of systematics
29  std::vector<const ISyst*> systs = getAllNusSysts();
30 
31  std::string folder = "/nova/ana/steriles/Ana01/";
32  std::string filenm = "FDDataMCSystBand";
33 
34  std::string loadLocation = folder + filenm + ".root";
35  std::string textLocation = folder + "RWithSysts.txt";
36 
37  TFile* rootL = new TFile(loadLocation.c_str(), "READ");
38  TDirectory* tmpL = gDirectory;
39  TDirectory* loadDir = gDirectory;
40 
41  loadDir->cd((loadLocation + ":/pred").c_str());
43 
44  loadDir->cd((loadLocation + ":/sData").c_str());
45  Spectrum sData = *Spectrum::LoadFrom(gDirectory);
46 
47  tmpL->cd();
48  rootL->Close(); // Don't forget to close the file!
49 
50  double r = 1.;
51  std::vector<double> rUps, rDns;
52  std::vector<double> rDiffUps, rDiffDns;
53 
57  r = CalcR(sData, nomCC, nomNC);
58 
59  double upSq = 0., dnSq = 0.;
60 
61  for(const ISyst* syst: systs) {
62  SystShifts shifts;
63 
64  shifts.SetShift(syst, +1);
67  double rUp = CalcR(sData, upCC, upNC);
68 
69  shifts.SetShift(syst, -1);
72  double rDn = CalcR(sData, dnCC, dnNC);
73 
74  if(rUp > r && rDn > r) {
75  rUp = std::max(rUp, rDn);
76  rDn = r;
77  }
78  if(rUp < r && rDn < r) {
79  rDn = std::min(rUp, rDn);
80  rUp = r;
81  }
82  if(rUp < r && rDn > r) {
83  std::swap(rUp, rDn);
84  }
85 
86  double diffUp = CalcDiffR(r, rUp);
87  double diffDn = CalcDiffR(r, rDn);
88 
89  rUps.push_back(rUp);
90  rDns.push_back(rDn);
91  rDiffUps.push_back(diffUp);
92  rDiffDns.push_back(diffDn);
93 
94  upSq += diffUp*diffUp;
95  dnSq += diffDn*diffDn;
96  }
97 
98  FILE* textF;
99  textF = fopen(textLocation.c_str(), "w");
100 
101  fprintf(textF, "R = %4.3f +%.3f -%.3f\n\n",
102  r, sqrt(upSq), sqrt(dnSq));
103  fprintf(textF, "Systematic - R Up, R Dn\n\n");
104 
105  for(int i = 0, n = systs.size(); i < n; ++i) {
106  fprintf(textF, "%9.9s -- %4.3f, %4.3f\n",
107  systs[i]->ShortName().c_str(), rUps[i], rDns[i]);
108  fprintf(textF, "R = %4.3f +%.3f -%.3f\n\n",
109  r, rDiffUps[i], rDiffDns[i]);
110  }
111 
112  fclose(textF);
113 }
114 
115 //----------------------------------------------------------------------
116 double CalcR(const Spectrum& sData, const Spectrum& sMCCC, const Spectrum& sMCNC)
117 {
118  TH1* hData = sData.ToTH1(sData.POT());
119  TH1* hMCCC = sMCCC.ToTH1(kSecondAnaPOT);
120  TH1* hMCNC = sMCNC.ToTH1(kSecondAnaPOT);
121 
122  double nData = hData->Integral();
123  double nMCCC = hMCCC->Integral();
124  double nMCNC = hMCNC->Integral();
125  return ((nData - nMCCC - 14.32527) / nMCNC);
126 }
127 
128 //----------------------------------------------------------------------
129 double CalcDiffR(double r, double rdiff)
130 {
131  return std::abs(r - rdiff);
132 }
T max(const caf::Proxy< T > &a, T b)
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
static std::unique_ptr< PredictionCombinePeriods > LoadFrom(TDirectory *dir, const std::string &name)
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
T sqrt(T number)
Definition: d0nt_math.hpp:156
Adapt the PMNS_Sterile calculator to standard interface.
osc::OscCalcDumb calc
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
float abs(float number)
Definition: d0nt_math.hpp:39
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
void CalcRWithSysts()
fclose(fg1)
Charged-current interactions.
Definition: IPrediction.h:39
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
Spectrum PredictComponentSyst(osc::IOscCalc *calc, const SystShifts &syst, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
std::vector< const ISyst * > getAllNusSysts()
Get a vector of all the nus group systs.
Definition: NusSysts.cxx:202
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
double CalcR(const Spectrum &sData, const Spectrum &sMCCC, const Spectrum &sMCNC)
Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
double POT() const
Definition: Spectrum.h:227
double CalcDiffR(double r, double rdiff)
Neutral-current interactions.
Definition: IPrediction.h:40
TRandom3 r(0)
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Sum MC predictions from different periods scaled according to data POT targets.
T min(const caf::Proxy< T > &a, T b)
All neutrinos, any flavor.
Definition: IPrediction.h:26
const double kSecondAnaPOT
Definition: Exposures.h:87
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:81
enum BeamMode string