CalcRWithSystsNus17.C
Go to the documentation of this file.
1 // Macro to calculate FD R-value
2 // TODO: Make it print out maximal mixing numbers too
3 
7 #include "CAFAna/Core/Spectrum.h"
10 #include "CAFAna/Systs/Systs.h"
11 
12 #include "OscLib/OscCalcPMNSOpt.h"
13 #include "OscLib/OscCalcDumb.h"
14 
15 #include "TFile.h"
16 #include "TH1.h"
17 
18 #include <iostream>
19 #include <string>
20 #include <vector>
21 
22 using namespace ana;
23 
24 
26 
27 
28 double CalcR(const Spectrum& sData, const Spectrum& sMCCC, const Spectrum& sMCNC, bool isHigh);
29 double CalcDiffR(double r, double rdiff);
30 
32 
33 void CalcRWithSystsNus17(std::string outname, bool isHigh, bool isUpper)
34 {
35  TH1::AddDirectory(0);
36 
37  // Create vector of systematics
38  std::vector<const ISyst*> systs = getAllNusFDSysts();
39 
40  std::string folder = "/nova/ana/steriles/Nus17/";
41  std::string filenm = "predInterpSysts_nus17_v1";
42 
43  std::string loadLocation = folder + filenm + ".root";
44  std::string textLocation = "nus17_RWithSysts.txt";
45 
46  TFile* rootL = new TFile(loadLocation.c_str(), "READ");
47  TDirectory* tmpL = gDirectory;
48  TDirectory* loadDir = gDirectory;
49 
50  loadDir->cd((loadLocation + ":/nus17_predI").c_str());
51  std::unique_ptr<PredictionInterp> pred = PredictionInterp::LoadFrom(gDirectory);
52 
53  Spectrum* intime = LoadFromFile<Spectrum>("$NOVA_ANA/steriles/Nus17/fout_nus17_box_opening_restricted.root","spec_nus_vise_numi").release();
54 
55  TH1 *hSpectrum = intime->ToTH1(intime->POT());
56  Spectrum sData(hSpectrum, intime->POT(), intime->Livetime());
57  std::cout << "Data spectrum: " << sData.POT()
58  <<" POT and "
59  << sData.ToTH1(intime->POT())->GetSumOfWeights()
60  << " events" << std::endl;
61 
62  TFile *cos = new TFile("$NOVA_ANA/steriles/Nus17/nus17_cosmic_background.root","read");
63  TH1 *hcosmic = (TH1*)cos->Get("cosall");
65  std::cout << "Cosmic Data spectrum: " << cosmic.POT()
66  <<" POT and "
67  << cosmic.ToTH1(cosmic.POT())->GetSumOfWeights()
68  << " events"
69  << ", Livetime: "
70  << cosmic.Livetime() << std::endl;
71  delete hcosmic;
72  cos->Close();
73 
74 
75  tmpL->cd();
76  rootL->Close(); // Don't forget to close the file!
77 
78  double r = 1.;
79  std::vector<double> rUps, rDns;
80  std::vector<double> rDiffUps, rDiffDns;
81 
82  //osc::OscCalcSterile* calc = DefaultSterileCalc(3);
84  calc = MyDefaultOscCalc();
85  calc->SetTh23(asin(sqrt(0.404))); // 0.624, dcp= 0.74pi
86  calc->SetdCP(1.48*M_PI);
87  calc->SetDmsq32(2.67e-3);
88  if(isUpper) {
89  std::cout << "IS UPPER" << std::endl;
90  calc->SetTh23(asin(sqrt(0.624)));
91  calc->SetdCP(0.74*M_PI);
92  }
93 
94 
95  Spectrum nomNC = pred.get()->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth);
96  Spectrum nomCC = pred.get()->PredictComponent(calc, Flavors::kAll, Current::kCC, Sign::kBoth);
97  Spectrum nomCCnue = pred.get()->PredictComponent(calc, Flavors::kAllNuE, Current::kCC, Sign::kBoth);
98  Spectrum nomCCnumu = pred.get()->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth);
99  Spectrum nomCCnutau = pred.get()->PredictComponent(calc, Flavors::kAllNuTau, Current::kCC, Sign::kBoth);
100 
101  TH1* hMCNC = nomNC.ToTH1(kAna2017POT);
102  TH1* hMCCCnue = nomCCnue.ToTH1(kAna2017POT);
103  TH1* hMCCCnumu = nomCCnumu.ToTH1(kAna2017POT);
104  TH1* hMCCCnutau = nomCCnutau.ToTH1(kAna2017POT);
105  double nMCNC=0.;
106  double nMCCCnue=0.;
107  double nMCCCnumu=0.;
108  double nMCCCnutau=0.;
109  double nData = 0.;
110  double nCos = 0.;
111  nData = sData.ToTH1(sData.POT())->Integral(0,25);
112  nCos = cosmic.ToTH1(cosmic.POT())->Integral(0,25);
113  nMCNC = hMCNC->Integral(0, 25);
114  nMCCCnue = hMCCCnue->Integral(0, 25);
115  nMCCCnumu = hMCCCnumu->Integral(0, 25);
116  nMCCCnutau = hMCCCnutau->Integral(0, 25);
117  std::cout << "Data: " << nData << std::endl;
118  std::cout << "Cosmic: " << nCos << std::endl;
119  std::cout << "NC: " << nMCNC << std::endl;
120  std::cout << "CCnue: " << nMCCCnue << std::endl;
121  std::cout << "CCnumu: " << nMCCCnumu << std::endl;
122  std::cout << "CCnutau: " << nMCCCnutau << std::endl;
123 
124  r = CalcR(sData, nomCC, nomNC, isHigh);
125 
126  double upSq = 0., dnSq = 0.;
127 
128  for(const ISyst* syst: systs) {
129  SystShifts shifts;
130 
131  shifts.SetShift(syst, +1);
132  Spectrum upNC = pred.get()->PredictComponentSyst(calc, shifts, Flavors::kAll, Current::kNC, Sign::kBoth);
133  Spectrum upCC = pred.get()->PredictComponentSyst(calc, shifts, Flavors::kAll, Current::kCC, Sign::kBoth);
134  double rUp = CalcR(sData, upCC, upNC, isHigh);
135 
136  shifts.SetShift(syst, -1);
137  Spectrum dnNC = pred.get()->PredictComponentSyst(calc, shifts, Flavors::kAll, Current::kNC, Sign::kBoth);
138  Spectrum dnCC = pred.get()->PredictComponentSyst(calc, shifts, Flavors::kAll, Current::kCC, Sign::kBoth);
139  double rDn = CalcR(sData, dnCC, dnNC, isHigh);
140 
141  if(rUp > r && rDn > r) {
142  rUp = std::max(rUp, rDn);
143  rDn = r;
144  }
145  if(rUp < r && rDn < r) {
146  rDn = std::min(rUp, rDn);
147  rUp = r;
148  }
149  if(rUp < r && rDn > r) {
150  std::swap(rUp, rDn);
151  }
152 
153  double diffUp = CalcDiffR(r, rUp);
154  double diffDn = CalcDiffR(r, rDn);
155 
156  rUps.push_back(rUp);
157  rDns.push_back(rDn);
158  rDiffUps.push_back(diffUp);
159  rDiffDns.push_back(diffDn);
160 
161  upSq += diffUp*diffUp;
162  dnSq += diffDn*diffDn;
163  }
164 
165  FILE* textF;
166  textF = fopen(outname.c_str(), "w");
167 
168  fprintf(textF, "R = %4.3f +%.3f -%.3f\n\n",
169  r, sqrt(upSq), sqrt(dnSq));
170  fprintf(textF, "Systematic - R Up, R Dn\n\n");
171 
172  for(int i = 0, n = systs.size(); i < n; ++i) {
173  fprintf(textF, "%9.9s -- %4.3f, %4.3f\n",
174  systs[i]->ShortName().c_str(), rUps[i], rDns[i]);
175  fprintf(textF, "R = %4.3f +%.3f -%.3f\n\n",
176  r, rDiffUps[i], rDiffDns[i]);
177  }
178 
179  fclose(textF);
180 }
181 
182 //----------------------------------------------------------------------
183 double CalcR(const Spectrum& sData, const Spectrum& sMCCC, const Spectrum& sMCNC, bool isHigh)
184 {
185 
186  TH1* hData = sData.ToTH1(sData.POT());
187  TH1* hMCCC = sMCCC.ToTH1(kAna2017POT);
188  TH1* hMCNC = sMCNC.ToTH1(kAna2017POT);
189 
190  double nData = 0.;
191  double nMCCC = 0.;
192  double nMCNC = 0.;
193  double Rvalue = 0.;
194  if(isHigh){
195  nData = hData->Integral(25, 100);
196  nMCCC = hMCCC->Integral(25, 100);
197  nMCNC = hMCNC->Integral(25, 100);
198  Rvalue = ((nData - nMCCC - 7.89788) / nMCNC);
199  }
200  else{
201  nData = hData->Integral(0, 25);
202  nMCCC = hMCCC->Integral(0, 25);
203  nMCNC = hMCNC->Integral(0, 25);
204  Rvalue = ((nData - nMCCC) / nMCNC);
205 
206  }
207  return Rvalue;
208 }
209 
210 //----------------------------------------------------------------------
211 double CalcDiffR(double r, double rdiff)
212 {
213  return std::abs(r - rdiff);
214 }
215 
217 {
218  ResetOscCalcToDefault(calc);
219 }
220 
222 {
225  return ret;
226 }
T max(const caf::Proxy< T > &a, T b)
TCut intime("tave>=217.0 && tave<=227.8")
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
_OscCalcPMNSOpt< double > OscCalcPMNSOpt
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
void ResetOscCalcToMyDefault(osc::IOscCalcAdjustable *calc)
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
double Integral(const Spectrum &s, const double pot, int cvnregion)
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
T sqrt(T number)
Definition: d0nt_math.hpp:156
double CalcR(const Spectrum &sData, const Spectrum &sMCCC, const Spectrum &sMCNC, bool isHigh)
const double kAna2017POT
Definition: Exposures.h:174
osc::IOscCalcAdjustable * MyDefaultOscCalc()
osc::OscCalcDumb calc
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
virtual void SetDmsq32(const T &dmsq32)=0
#define M_PI
Definition: SbMath.h:34
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)
fclose(fg1)
Charged-current interactions.
Definition: IPrediction.h:39
const double kAna2017Livetime
Definition: Exposures.h:200
Sum up livetimes from individual cosmic triggers.
double CalcDiffR(double r, double rdiff)
double POT() const
Definition: Spectrum.h:227
void CalcRWithSystsNus17(std::string outname, bool isHigh, bool isUpper)
OStream cout
Definition: OStream.cxx:6
static std::unique_ptr< PredictionInterp > LoadFrom(TDirectory *dir, const std::string &name)
virtual void SetTh23(const T &th23)=0
Neutral-current interactions.
Definition: IPrediction.h:40
T cos(T number)
Definition: d0nt_math.hpp:78
TRandom3 r(0)
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
T min(const caf::Proxy< T > &a, T b)
All neutrinos, any flavor.
Definition: IPrediction.h:26
Float_t e
Definition: plot.C:35
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:230
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:81
virtual void SetdCP(const T &dCP)=0
T asin(T number)
Definition: d0nt_math.hpp:60
enum BeamMode string