sensitivity_tools.h
Go to the documentation of this file.
1 #pragma once
2 
14 #include "CAFAna/Systs/XSecSysts.h"
15 #include "CAFAna/Vars/FitVars.h"
16 
17 #include "CAFAna/Analysis/Style.h"
19 
20 #include "TCanvas.h"
21 #include "TBox.h"
22 #include "TLatex.h"
23 #include "TColor.h"
24 #include "TGraph.h"
25 #include "TVectorD.h"
26 #include "TF1.h"
27 #include "TLegend.h"
28 #include "TLine.h"
29 #include "TMarker.h"
30 #include "TStyle.h"
31 #include "TSystem.h"
32 #include "TGaxis.h"
33 
34 using namespace ana;
35 
36 #include "OscLib/IOscCalc.h"
37 
38 #include "TFile.h"
39 
42  bool corrSysts,
44  bool isFakeData,
45  bool minimizeMemory = false){
46 
48  if(decomp == "noextrap") extrap = kNoExtrap;
49  if(decomp == "prop") extrap = kProportional;
50  if(decomp == "combo") extrap = kCombo;
51  if(decomp == "fakexp") extrap = kFake;//Hack
52 
53  PredictionSystJoint2018 * pred_fid;
54 
55  std::string nue_pred_path = "";
56 
57  if(corrSysts){
58  std::cout << "Getting Nue Prediction from official predictions" << endl;
59  if(beam=="rhc") pred_fid = new PredictionSystJoint2018 (extrap, calc,
61  beam, isFakeData, true,
62  nue_pred_path, minimizeMemory);
63  else if(beam=="fhc") pred_fid = new PredictionSystJoint2018 (extrap, calc,
65  beam, isFakeData, true,
66  nue_pred_path, minimizeMemory);
67  else{ std::cout << "Beam mode needs to be either fhc or rhc" << endl; exit(1);}
68  }
69  else pred_fid = new PredictionSystJoint2018 (extrap, calc,{}, beam, isFakeData,
70  true, nue_pred_path, minimizeMemory);
71 
72 
73 
74  double pot;
75  if(beam=="rhc") pot = kFutureRHCPOT;
76  else pot = kFutureFHCPOT;
77  std::cout << "Scale to "<<pot<<std::endl;
78 
79  std::cout << "Nue prediction (no rock)"
80  << pred_fid->Predict(calc).Integral(pot) <<std::endl;
81  // Load rock
82 
83  double nonsScale = 1./10.45;
84  double swapScale = 1./12.91;
85 
86  std::string dirName = "/nova/ana/nu_e_ana/Ana2018/Predictions/FDRock/";
87  std::string fName2 = "fdrock_nue2018.root";
88  std::string rName;
89 
90  if(beam=="rhc") rName = "pred_FDRock_merg_rhc";
91  else rName = "pred_FDRock_merg_fhc";
92  if( decomp == "noextrap"){
93  if(beam=="rhc") rName = "pred_FDRock_rhc";
94  else rName = "pred_FDRock_fhc";
95  }
96 
97  auto rock = LoadFromFile<PredictionNoExtrap>(dirName + fName2, rName).release();
98 
99  std::cout << "Nue prediction (rock only - not scaled) "
100  << rock->Predict(calc).Integral(pot) <<std::endl;
101  std::cout << "Adding Rock" << endl;
102  auto nue_pred = new PredictionAddRock (pred_fid, rock, nonsScale, swapScale);
103 
104  std::cout << "Nue prediction (rock added) "
105  << nue_pred->Predict(calc).Integral(pot) << std::endl;
106  std::cout << "Pure rock events: "
107  << nue_pred->Predict(calc).Integral(pot) - pred_fid->Predict(calc).Integral(pot) << std::endl;
108 
109 
110  return nue_pred;
111 }
112 
113 //////////////////////////////////////////////////////////////////////
114 
115 std::pair <TH1D*, double > GetNueCosmics2018 (std::string beam, bool IsFuture = false){
116  Spectrum* outtime;
117 
119  std::string anaDir = "/nova/ana/nu_e_ana/Ana2018/";
120  if(beam=="rhc") name = "cosmic_spect_rhc";
121  else name="cosmic_spect_fhc";
122  outtime = LoadFromFile<Spectrum>
123  (anaDir+"/Predictions/cosmic/cosmic_prediction_real_data.root", name ).release();
124 
125  double livetime;
126  if(beam=="rhc") livetime = IsFuture? kFutureRHCLivetime: kFutureRHCLivetime ;
127  else livetime = IsFuture? kFutureFHCLivetime: kFutureFHCLivetime ;
128 
129  double outN = outtime->Integral(livetime, 0, kLivetime);
130 
131  std::cout << "\nAdding " << outtime->Integral(livetime, 0, kLivetime)
132  << " cosmics from out-of-time." << std::endl;
133  std::cout << "Which scale to " << outtime->Integral(livetime, 0,
134  kLivetime)
135  << " in the spill window." << std::endl;
136 
137  std::cout << "uncertainty " << 1/sqrt(outN) << std::endl;
138 
139  TH1 *h = outtime->ToTH1(livetime, kBlack, kSolid, kLivetime);
140 
141  return {outtime->ToTH1(livetime, kBlack, kSolid, kLivetime), 1/sqrt(outN)};
142 }
143 
144 //////////////////////////////////////////////////////////////////////
145 
147  Spectrum* intime;
148 
149  if(beam=="fhc"){
150  intime = LoadFromFile<Spectrum>
151  ("/nova/ana/nu_e_ana/Ana2018/Data2018/spectra.root", "data_spect_fhc").release();
152  }
153  else if(beam=="rhc"){
154  intime = LoadFromFile<Spectrum>
155  ("/nova/ana/nu_e_ana/Ana2018/Data2018/spectra.root", "data_spect_rhc").release();
156  }
157  else{ std::cout << "Beam mode != fhc or rhc" << std::endl; exit(1);}
158  double POT, livetime;
159  if(beam == "rhc") {livetime = kFutureRHCLivetime; POT = kFutureRHCPOT;}
160  if(beam == "fhc") {livetime = kFutureFHCLivetime; POT = kFutureFHCPOT;}
161 
162  std::cout << "Loaded data spectrum:\n"
163  << "POT " << intime->POT()
164  << " (expected " << POT << ")\n"
165  << "Livetime " << intime->Livetime()
166  << " (expected " << livetime << ")\n"
167  << "Integral " << intime->Integral(intime->POT())
168  << std::endl;
169 
170  return intime;
171 }
172 
173 //////////////////////////////////////////////////////////////////////
174 
175 std::vector <const IPrediction*> GetNumuPredictions2018(const int nq = 4,
176  bool useSysts = true,
177  std::string beam="fhc",
178  ENu2018ExtrapType numuExtrap = kNuMu,
179  bool minimizeMemory = false)
180 {
181  auto calc = DefaultOscCalc();
182 
183  std::vector <const IPrediction*> numu_preds;
184 
186  std::string anaDir = "/nova/ana/nu_mu_ana/Ana2018/";
187  if(beam=="rhc") filename=anaDir + "/Predictions/pred_numuconcat_rhc__numu2018.root";
188  else filename= anaDir + "/Predictions/pred_numuconcat_fhc__numu2018.root";
189 
190  std::cout << "\nLoading Numu Predictions\n" ;
191  for (int quant = 0; quant < nq; ++quant){
192  if(!useSysts) {
193  numu_preds.push_back(new PredictionSystJoint2018(numuExtrap, calc, beam, quant+1,
194  {}, filename, minimizeMemory));
195  }
196  else {
197  numu_preds.push_back(new PredictionSystJoint2018(numuExtrap, calc, beam, quant+1,
199  filename, minimizeMemory));
200  }
201  }
202  return numu_preds;
203 }
204 //////////////////////////////////////////////////////////////////////
205 
206 std::vector <std::pair <TH1D*, double> > GetNumuCosmics2018(const int nq = 4, std::string beam="fhc", double fromLivetime = kAna2018FHCLivetime, double toLivetime = kFutureFHCLivetime)
207 {
208  cout<<"Beam is "<<beam<<endl;
209  //cosmics
210  //the versionn in cafana has float not double and it's a pain
211  TFile *fcosm;
212 
213  std::string dir = "/nova/ana/nu_mu_ana/Ana2018/Cosmics/";
215  if (beam=="rhc") filename="cosmics_rhc__numu2018.root";
216  else filename="cosmics_fhc__numu2018.root";
217  fcosm = TFile::Open((dir+filename).c_str());
218 
219 
220  if(fcosm->IsZombie()) {std::cerr<< "bad cosmics\n"; exit(1);}
221  std::vector <std::pair <TH1D*, double > > numu_cosmics;
222 
223  for (int i = 1; i <=nq; ++i) {
224  auto h = (TH1D*) fcosm->Get((TString)"cosmics_q"+std::to_string(i))
225  ->Clone(UniqueName().c_str());
226 
227  h->Scale(toLivetime/fromLivetime);
228 // Spectrum* tmp = new Spectrum(h, 0, toLivetime);
229  numu_cosmics.push_back({h,1/sqrt(h->Integral())});
230  std::cout << "Cosmics all periods quantile " << i << " "
231  << h->Integral()
232  << " scale error "
233  //<< hs->GetBinContent(1)
234  << 1/sqrt(h->Integral())
235  << "\n";
236  }//end quantiles
237  return numu_cosmics;
238 }
239 //////////////////////////////////////////////////////////////////////
240 
241 std::vector <Spectrum * > GetNumuData2018(const int nq = 4, std::string beam = "fhc"){
242 
244  if(beam == "fhc") {file="/nova/ana/nu_mu_ana/Ana2018/Data/fd_dataspectrum_fhc_full__numu2018.root";}
245  if(beam == "rhc") {file="/nova/ana/nu_mu_ana/Ana2018/Data/fd_dataspectrum_rhc_full__numu2018.root";}
246 
247  std::vector <Spectrum *> numu_data;
248 
249  for (int i = 1; i <=nq; ++i) {
250 
251  auto f= new TFile(file.c_str());
252  if(f->IsZombie()) {std::cerr<< "bad data\n"; exit(1);}
253  numu_data.push_back(Spectrum::LoadFrom(f, ("fd_data_q"+std::to_string(i)).c_str())).release();
254 
255  }//end quantiles
256  return numu_data;
257 }
258 
259 //////////////////////////////////////////////////////////////////////
260 
262  const double pot, TH1D* cosmics, const double livetime)
263 {
264  auto temp = pred->Predict(calc).FakeData(pot);
265  if(cosmics)
266  temp += Spectrum(cosmics, pot, livetime);
267  return new Spectrum(temp);
268 }
269 
270 //////////////////////////////////////////////////////////////////////
271 
273  double dmsq32 = -5, double dCP_Pi = -5)
274 {
275  if(dCP_Pi > -1) kFitDeltaInPiUnits.SetValue(calc, dCP_Pi);
276  if(ssth23 > -1) kFitSinSqTheta23.SetValue(calc, ssth23);
277  if(dmsq32 > -1) kFitDmSq32.SetValue(calc, dmsq32);
278 }
279 
280 //////////////////////////////////////////////////////////////////////
281 
282 std::vector<const ISyst*> GetJointFitSystematicList(bool corrSysts, bool nueExclusive = false, bool numuExclusive=false, bool isFHC = true, bool isRHC = true , bool intersection = true){
283  std::vector<const ISyst*> ret ={} ;
284  if(corrSysts){
285  if(! intersection){
286  if(nueExclusive){
287 
288  ret.push_back(&kRadCorrNue);
289  ret.push_back(&kRadCorrNuebar);
290  ret.push_back(&k2ndClassCurrs);
293 
294  }
295  if(numuExclusive){
296 
300  if(isFHC) ret.push_back(&kAna2018NormRHC);
301  if(isRHC) ret.push_back(&kAna2018NormFHC);
303 
304  }
305  }
306  else{
307  if(nueExclusive) ret = getAllAna2018Systs(kNueAna2018);
308  if(numuExclusive) ret = getAllAna2018Systs(kNumuAna2018);
309  if(!nueExclusive && !numuExclusive) ret = getAllAna2018Systs(kJointAna2018, true, kBoth, true);
310  }
311  }
312  std::cout << "Return list of systematics: " ;
313  for(auto & s:ret) std::cout << s->ShortName() << ", ";
314  std::cout << std::endl;
315  return ret;
316 }
317  std::vector<std::pair<const ISyst*,const ISyst*> > GetCorrelations (bool isNue, bool isFHC){
318  std::vector<std::pair<const ISyst*,const ISyst*> > ret ={};
319  std::cout<<(isNue?"not nue ":"not numu ")<<(isFHC&&!isNue?" ":(isFHC?" FHC ":" RHC "));
320  auto badsyst = GetJointFitSystematicList(true, !isNue, isNue, isFHC, !isFHC, false);
321  for (auto & s:badsyst) ret.push_back ({s,NULL});
322 
323  return ret;
324 }
325 
326 /////////////////////////////////////////////////////////////////////
327 
328 // For slices plots, plot by true octant, and make separate plots of each
330 {
331 public:
332  virtual double GetValue(const osc::IOscCalcAdjustable* osc) const {
333  return util::sqr(sin(osc->GetTh23()));};
334  virtual void SetValue(osc::IOscCalcAdjustable* osc, double val) const {
335  osc->SetTh23(asin(sqrt(Clamp(val))));};
336  virtual std::string ShortName() const {return "ssth23LO";}
337  virtual std::string LatexName() const {return "sin^{2}#theta_{23} LO";}
338  virtual double LowLimit() const {return 0.01;}
339  virtual double HighLimit() const {return 0.499;}
340 };
341 
343 
345 {
346 public:
347  virtual double GetValue(const osc::IOscCalcAdjustable* osc) const {
348  return util::sqr(sin(osc->GetTh23()));};
349  virtual void SetValue(osc::IOscCalcAdjustable* osc, double val) const {
350  osc->SetTh23(asin(sqrt(Clamp(val))));};
351  virtual std::string ShortName() const {return "ssth23UO";}
352  virtual std::string LatexName() const {return "sin^{2}#theta_{23} UO";}
353  virtual double LowLimit() const {return 0.501;}
354  virtual double HighLimit() const {return 0.99;}
355 };
357 
359 {
360 public:
361  virtual double GetValue(const osc::IOscCalcAdjustable* osc) const {
362  return util::sqr(sin(osc->GetTh23()));};
363  virtual void SetValue(osc::IOscCalcAdjustable* osc, double val) const {
364  osc->SetTh23(asin(sqrt(Clamp(val))));};
365  virtual std::string ShortName() const {return "ssth23MaxMix";}
366  virtual std::string LatexName() const {return "sin^{2}#theta_{23} MaxMix";}
367  virtual double LowLimit() const {return 0.49;}
368  virtual double HighLimit() const {return 0.51;}
369 };
371 
372 ////////////////////////////////////////////////////////////////////////////
373 //////////////////////// Slice, Contours Style /////////////////////////////
374 ///////////////////////////////////////////////////////////////////////////
375 //move functions from plotting script here
376 void latex(double x, double y, std::string txt, double ang = 0, int align = 22, double size = 2/30.)
377 {
378  TLatex* ltx = new TLatex(x, y, txt.c_str());
379  ltx->SetNDC();
380  ltx->SetTextSize(size);
381  ltx->SetTextAlign(align);
382  ltx->SetTextAngle(ang);
383 cout<<"i can draw!"<<endl;
384  ltx->Draw();
385 }
386 
387 void DrawSliceCanvas(bool IsTh23, TString slicename, const double ymax,
388  const double xmin = 0, const double xmax = 2.){
389  auto c1 = new TCanvas(ana::UniqueName().c_str());
390  c1->SetFillStyle(0);
391  c1->SetLeftMargin(0.14);
392  c1->SetBottomMargin(0.15);
393  TString title;
394  TString xytitles = IsTh23?"sin^{2}#theta_{23};":"#delta_{CP} / #pi;";//Significance (#sqrt{#Delta#chi^{2}})";
395 
396  if(slicename.Contains("mh")) {
397  title = ";" + xytitles + "Significance of hierarchy resolution(#sigma)";
398 // title = "Rejection of mass hierarchy;" + xytitles;
399  }
400  if(slicename.Contains("oct")) {
401  title = ";" + xytitles + "Significance of octant determination(#sigma)";
402  }
403  if(slicename.Contains("maxmix")){
404  title = ";" + xytitles + "Significance of maximal mixing rejection(#sigma)";
405  gStyle->SetTitleSize(0.045, "y");
406  }
407  if(slicename.Contains("dcp")){
408  title = ";" + xytitles + "Significance of CP violation(#sigma)";
409  }
410  if(slicename.Contains("dmsq")) {
411  title = ";|#Deltam^{2}_{32}| (10^{-3} eV^{2});Significance (#sqrt{#Delta#chi^{2}})";
412  }
413  if(slicename.Contains("th13")) {
414  title = ";sin^{2}2#theta_{13};Significance (#sqrt{#Delta#chi^{2}})";
415  }
416  auto axes = new TH2F("",title, 100, xmin, xmax, 100, 0, ymax);
417  if(slicename.Contains("fccorr")) axes->GetYaxis()->SetTitle("Significance (#sigma)");
418  axes->GetXaxis()->CenterTitle();
419  axes->GetYaxis()->CenterTitle();
420  axes->GetXaxis()->SetTitleOffset(0.9);
421  axes->Draw();
422 // if(!IsTh23) XAxisDeltaCPLabels(axes);
423  c1->RedrawAxis();
424 
425  }
426 
427 TLegend * SliceLegend(std::vector <std::pair <TGraph*, TString > > & graphs, bool isDelta, double low=0.68, bool syst=false)
428 {
429  TLegend * leg;
430  leg = new TLegend(0.42, low, 0.58, 0.85);
431 
432  leg->SetTextSize(1.55/30.);
433  leg->SetFillStyle(0);
434  for (auto &gr:graphs){
435  leg->AddEntry(gr.first,gr.second, "l");
436  }
437  leg->SetHeader("Stat. Only");
438  if(syst) leg->SetHeader("2018 Syst.");
439  return leg;
440 }
441 
442 void DrawContourCanvas (TString surfName, double minx = 0,
443  double maxx = 2, double miny = 0, double maxy = 1){
444  auto c1 = new TCanvas(ana::UniqueName().c_str());
445  c1->SetFillStyle(0);
446  c1->SetLeftMargin(0.14);
447  c1->SetBottomMargin(0.15);
448  TH2* axes = new TH2F();
449  TString title;
450  if(surfName.Contains("delta_th23")) title=";#delta_{CP};sin^{2}#theta_{23}";
451  if(surfName.Contains("th13_delta")) title=";sin^{2}2#theta_{13};#delta_{CP}/#pi";
452  if(surfName.Contains("th23_dm32")) title=";sin^{2}#theta_{23};#Deltam^{2}_{32} (10^{-3}eV^{2})";
453  axes = new TH2F("",title,100,minx,maxx,100,miny,maxy);
454 // CenterTitles(axes);
455  axes->GetXaxis()->CenterTitle();
456  axes->GetYaxis()->CenterTitle();
457  axes->Draw();
458  axes->GetXaxis()->SetTitleSize(kBlessedTitleFontSize);
459  axes->GetYaxis()->SetTitleSize(kBlessedTitleFontSize);
460  axes->GetXaxis()->SetLabelSize(kBlessedLabelFontSize);
461  axes->GetYaxis()->SetLabelSize(kBlessedLabelFontSize);
462  axes->GetXaxis()->SetTitleOffset(0.8);
463  axes->GetYaxis()->SetTitleOffset(0.8);
464  TGaxis::SetMaxDigits(3);
465  if(surfName.Contains("th23_dm32")) {
466  axes->GetYaxis()->SetDecimals() ;
467  axes->GetYaxis()->SetTitleOffset(0.85);
468  axes->GetYaxis()->SetLabelOffset(0.002);
469  axes->GetYaxis()->SetLabelSize(0.058);
470  axes->GetYaxis()->SetTitleSize(0.078);
471  }
472  if(surfName.Contains("delta_th23")) XAxisDeltaCPLabels(axes);
473  c1->RedrawAxis();
474 }
475 
476 
477 
Spectrum * GetNueData2018(std::string beam)
const XML_Char * name
Definition: expat.h:151
Simple analyzer module to print to text file information required for c++ based block alignment study...
TCut intime("tave>=217.0 && tave<=227.8")
const NOvARwgtSyst k2ndClassCurrs("2ndclasscurr","Second class currents", novarwgt::kSimpleSecondClassCurrentsSystKnob)
Second-class current syst. See documentation in NOvARwgt.
Definition: XSecSysts.h:71
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::map< std::string, double > xmax
double ssth23
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
Loads shifted spectra from files.
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:18
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
double maxy
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Definition: FitVars.cxx:177
void XAxisDeltaCPLabels(TH1 *axes, bool t2kunits)
Label the x-axis with fractions of pi.
T sqrt(T number)
Definition: d0nt_math.hpp:156
const FitSinSqTheta23UpperOctant kFitSinSqTheta23UpperOctant
Definition: FitVars.cxx:16
General interface to oscillation calculators.
Definition: StanTypedefs.h:22
const IConstrainedFitVar * kFitSinSqTheta23MaxMix
virtual double HighLimit() const
OStream cerr
Definition: OStream.cxx:7
virtual double LowLimit() const
const double kFutureRHCPOT
Definition: Exposures.h:245
virtual T GetTh23() const
Definition: IOscCalc.h:94
virtual double GetValue(const osc::IOscCalcAdjustable *osc) const
string filename
Definition: shutoffs.py:106
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Forward to wrapped Var&#39;s SetValue()
void AddJointAna2018NotCorrelSysts(std::vector< const ISyst * > &systs, const EAnaType2018 ana, const BeamType2018 beam)
const double kFutureFHCPOT
Definition: Exposures.h:244
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:249
virtual double HighLimit() const
virtual double LowLimit() const
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const NOvARwgtSyst kRadCorrNuebar("radcorrnuebar","Radiative corrections for #bar{#nu}_{e}", novarwgt::kSimpleRadiativeCorrNuebarXsecSystKnob)
Radiative corrections syst (nuebars). See documentation in NOvARwgt.
Definition: XSecSysts.h:67
Double_t ymax
Definition: plot.C:25
Spectrum * GetFakeData(const IPrediction *pred, osc::IOscCalc *calc, const double pot, TH1D *cosmics, const double livetime)
const XML_Char * s
Definition: expat.h:262
virtual double GetValue(const osc::IOscCalcAdjustable *osc) const
std::vector< const ISyst * > getAllAna2018Systs(const EAnaType2018 ana, const bool smallgenie, const BeamType2018 beam, bool isFit)
std::pair< TH1D *, double > GetNueCosmics2018(std::string beam, bool IsFuture=false)
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
std::vector< const IPrediction * > GetNumuPredictions2018(const int nq=4, bool useSysts=true, std::string beam="fhc", ENu2018ExtrapType numuExtrap=kNuMu, bool minimizeMemory=false)
Spectrum FakeData(double pot) const
Synonymous with AsimovData(). Retained for compatibility.
Definition: Spectrum.cxx:349
virtual void SetValue(osc::IOscCalcAdjustable *osc, double val) const
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
const double kFutureFHCLivetime
Definition: Exposures.h:247
#define pot
void SetFakeCalc(osc::IOscCalcAdjustable *calc, double ssth23=-5, double dmsq32=-5, double dCP_Pi=-5)
TLegend * SliceLegend(std::vector< std::pair< TGraph *, TString > > &graphs, bool isDelta, double low=0.68, bool syst=false)
T Clamp(T val) const
Definition: IFitVar.h:60
void DrawSliceCanvas(bool IsTh23, TString slicename, const double ymax, const double xmin=0, const double xmax=2.)
std::vector< float > Spectrum
Definition: Constants.h:610
const double kFutureRHCLivetime
Definition: Exposures.h:248
double POT() const
Definition: Spectrum.h:227
virtual std::string LatexName() const
static bool isFHC
Oscillation probability calculators.
Definition: Calcs.h:5
virtual std::string ShortName() const
void DrawContourCanvas(TString surfName, double minx=0, double maxx=2, double miny=0, double maxy=1)
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
virtual std::string LatexName() const
const IPrediction * GetNuePrediction2018(std::string decomp, osc::IOscCalc *calc, bool corrSysts, std::string beam, bool isFakeData, bool minimizeMemory=false)
void latex(double x, double y, std::string txt, double ang=0, int align=22, double size=2/30.)
Spectrum Predict(osc::IOscCalc *calc) const override
T sin(T number)
Definition: d0nt_math.hpp:132
const DummyAnaSyst kAna2018NormRHC("NormRHC2018","RHC. Norm.")
std::string dirName
Definition: PlotSpectra.h:47
virtual double LowLimit() const
double dmsq32
TDirectory * dir
Definition: macro.C:5
double livetime
Definition: saveFDMCHists.C:21
const FitSinSqTheta23LowerOctant kFitSinSqTheta23LowerOctant
Definition: FitVars.cxx:17
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
virtual void SetTh23(const T &th23)=0
void AddNonLoadable2018Systs(std::vector< const ISyst * > &systs, const EAnaType2018 ana)
exit(0)
virtual double HighLimit() const
TCut outtime("(tave>0.0&&tave<200.0) || (tave>250.0&&tave<450.0)")
TFile * file
Definition: cellShifts.C:17
std::vector< std::pair< TH1D *, double > > GetNumuCosmics2018(const int nq=4, std::string beam="fhc", double fromLivetime=kAna2018FHCLivetime, double toLivetime=kFutureFHCLivetime)
c1
Definition: demo5.py:24
std::vector< Spectrum * > GetNumuData2018(const int nq=4, std::string beam="fhc")
std::vector< std::pair< const ISyst *, const ISyst * > > GetCorrelations(bool isNue, bool isFHC, bool ptExtrap)
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const Float_t kBlessedLabelFontSize
Definition: Style.h:90
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
virtual std::string ShortName() const
const Float_t kBlessedTitleFontSize
Definition: Style.h:89
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
std::vector< const ISyst * > GetJointFitSystematicList(bool corrSysts, bool nueExclusive=false, bool numuExclusive=false, bool isFHC=true, bool isRHC=true, bool intersection=true, bool ptExtrap=true)
const NOvARwgtSyst kRadCorrNue("radcorrnue","Radiative corrections for #nu_{e}", novarwgt::kSimpleRadiativeCorrNueXsecSystKnob)
Radiative corrections syst (nues). See documentation in NOvARwgt.
Definition: XSecSysts.h:64
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:230
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
void rock(std::string suffix="full")
Definition: rock.C:28
virtual void SetValue(osc::IOscCalcAdjustable *osc, double val) const
const DummyAnaSyst kAna2018NormFHC("NormFHC2018","FHC. Norm.")
const double kAna2018FHCLivetime
Definition: Exposures.h:213
virtual std::string LatexName() const
virtual double GetValue(const osc::IOscCalcAdjustable *osc) const
virtual void SetValue(osc::IOscCalcAdjustable *osc, double val) const
T asin(T number)
Definition: d0nt_math.hpp:60
virtual std::string ShortName() const
enum BeamMode string