plot_nue_xsec_pred.C
Go to the documentation of this file.
2 #include "CAFAna/Core/ISyst.h"
3 #include "CAFAna/Systs/XSecSystLists.h" // getAllXsecSysts_2017
7 #include "CAFAna/Core/Spectrum.h"
10 #include "CAFAna/Analysis/Style.h"
11 #include "CAFAna/Analysis/Plots.h"
12 //#include "CAFAna/nue/SecondAna/draw_plots_util.h"
13 // TODO: #include rootlogon
14 
15 #include "OscLib/IOscCalc.h"
16 
18 
20 
21 #include "TString.h"
22 
23 using namespace ana;
24 
25 void plot_nue_xsec_pred(TString suffix="test", int startIdx=0, int stopIdx=-1)
26 {
27  // Set Oscillators
30 
31  calc->SetDmsq32(2.44e-3);
32  calc->SetTh13(asin(sqrt(.085))/2);
33  calc->SetTh23(asin(sqrt(0.45)));
34  calc->SetdCP(1.59*TMath::Pi());
35 
36  // Set POT
37  double pot = 9.e20;
38 
39  // Filename for output prediction file
40  TString filename = "pred_xsec_fhc_" + suffix + ".root";
41 
42  PredictionWriter* reader = new PredictionWriter(filename, "READ");
43 
44  // Set analysis name
45  TString analysisName = "pred_propDecomp";
46 
47  // Get vector of xsec systematics
48  std::vector<const ISyst*> systs = getAllXsecSysts_2017();
49  std::cout << "[make_nue_xsec_pred] Found " << systs.size()
50  << " xsec systematics!" << std::endl;
51  if (stopIdx < 0) // assume all systematics, set as length of vector
52  stopIdx = systs.size();
53  std::cout << "[make_nue_xsec_pred] Processing systematics:" << std::endl;
54  for (std::vector<const ISyst*>::iterator syst = systs.begin()+startIdx;
55  (syst != systs.begin()+stopIdx) && (syst != systs.end()); ++syst)
56  std::cout << std::distance(systs.begin(), syst) << ": "
57  << (*syst)->ShortName() << std::endl;
58 
59  // Set energyVar
60  TString energyVar = "EnergyCVN2D";
61 
62  // Just plot extrap and ppfx_xsec weights
63 
64  // Nominal
65  const TString predNomDir = reader->ReadPrediction(
66  analysisName, "Nominal", 0, true, energyVar, "ppfx_xsec");
67  PredictionExtrap* predNom = LoadFromFile<PredictionExtrap>(
68  reader->GetFilename().Data(), predNomDir.Data()).release();
69  Spectrum specNom = predNom->Predict(calc);
70  specNom.ToTH1(pot)->Draw();
71 
72  // Just a single plot for now, will add more now that it is working
73  // Loop over systematics
74  // Loop over systematics
75  for (std::vector<const ISyst*>::iterator syst = systs.begin()+startIdx;
76  (syst != systs.begin()+stopIdx) && (syst != systs.end()); ++syst)
77  {
78  const TString predPlusDir = reader->ReadPrediction(
79  analysisName, (*syst)->ShortName(), 2, true, energyVar, "ppfx_xsec");
80  PredictionExtrap* predPlus = LoadFromFile<PredictionExtrap>(
81  reader->GetFilename().Data(), predPlusDir.Data()).release();
82 
83  const TString predMinusDir = reader->ReadPrediction(
84  analysisName, (*syst)->ShortName(), -2, true, energyVar, "ppfx_xsec");
85  PredictionExtrap* predMinus = LoadFromFile<PredictionExtrap>(
86  reader->GetFilename().Data(), predMinusDir.Data()).release();
87 
88  std::vector<Spectrum> upShifts;
89  upShifts.emplace_back(predPlus->Predict(calc));
90  std::vector<Spectrum> downShifts;
91  downShifts.emplace_back(predMinus->Predict(calc));
92  //TGraphAsymmErrors* graph = PlotWithSystErrorBand(
93  // specNom, upShifts, downShifts, pot);
94 
95  //TCanvas* canvas = new TCanvas(systName);
96  //canvas->cd();
97  //graph->Draw();
98 
99  //TLegend* legend = new TLegend();
100  //legend->AddEntry(graph);
101  //legend->Draw();
102 
103  ComparePredictions(predNom, predPlus, predMinus,
104  analysisName.Data(), (*syst)->ShortName(), energyVar.Data(), 4);
105 
106  PredRatioToNom(predNom, predPlus, predMinus,
107  analysisName.Data(), (*syst)->ShortName(), energyVar.Data(), 4);
108  }
109 };
Pass neutrinos through unchanged.
Definition: IOscCalc.h:45
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
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
virtual void SetTh13(const T &th13)=0
const TString GetFilename() const
string filename
Definition: shutoffs.py:106
void plot_nue_xsec_pred(TString suffix="test", int startIdx=0, int stopIdx=-1)
unsigned distance(const T &t1, const T &t2)
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
virtual void SetDmsq32(const T &dmsq32)=0
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
void PredRatioToNom(PredictionExtrap *predPropMC, PredictionExtrap *predPropMCPlus, PredictionExtrap *predPropMCMinus, const std::string &predName, const std::string &systName, const std::string &varName, int varIdx)
Definition: NueSystFuncs.h:124
Spectrum Predict(osc::IOscCalc *calc) const override
#define pot
osc::NoOscillations noosc
Definition: NueSystFuncs.h:44
OStream cout
Definition: OStream.cxx:6
virtual void SetTh23(const T &th23)=0
const TString ReadPrediction(const TString &analysisName, const TString &systName, const int &sigma, const bool &extrap=true, const TString &energyVar="EnergyCVN2D", const TString &weight="kXSecCVWgt2017")
Take the output of an extrapolation and oscillate it as required.
Float_t e
Definition: plot.C:35
void ComparePredictions(IPrediction *pred_no, IPrediction *pred_sh, TString plottitle, TString out_name, TString pidtag, bool printtable=true, bool pidaxis=false)
std::vector< const ISyst * > getAllXsecSysts_2017()
Get master XSec syst list for 2017 analyses.
virtual void SetdCP(const T &dCP)=0
T asin(T number)
Definition: d0nt_math.hpp:60