make_prediction_extrap.C
Go to the documentation of this file.
1 #pragma once
2 
5 
6 using namespace ana;
7 
8 void make_prediction_extrap(std::string polarity = "rhc",
9  std::string period = "full",
10  bool nh = true,
11  bool systematics = false,
12  bool calc3a = true)
13 {
14 
15 
16  std::cout << "\n=================================" << std::endl;
17  std::cout << " make predictions extrap " << std::endl;
18  std::cout << "=================================\n" << std::endl;
19 
20  std::string extrapName = "extrap";
21 
22  // production files
23  std::string production = "";
24  if(polarity == "fhc"){
25  production += ".d";
26  }
27  if(polarity == "rhc"){
28  production += ".e";
29  }
30 
31  // systematics
32  std::string systName = "";
34  if(systematics){
35  systName = "systs";
36  }
37  if(!systematics){
38  systName = "stats";
39  systs = {};
40  }
41  for (auto & sys:systs) std::cout << sys->ShortName() << std::endl;
42 
43  // calculator
44  std::string calcName = "2018calc";
45  std::string hierarchy = "";
47  if(nh){
48  hierarchy = "nh";
49  calc->SetTh23(asin(sqrt(0.559455)));
50  calc->SetDmsq32(0.00244);
51  }
52  if(!nh){
53  hierarchy = "ih";
54  calc->SetTh23(asin(sqrt(0.559455)));
55  calc->SetDmsq32(-0.00244);
56  }
57 
58 
59  // selection and weight
60  const Cut kNumuNDCut = kNumuCutND2018;
61  const Cut kNumuFDCut = kNumuCutFD2018;
62  const Var kXSecWeight = kXSecCVWgt2018;
63  std::string OutDir = "";
64  std::cout << "polarity: " << polarity << ", period: " << period << ", files: " << production << std::endl;
65  std::cout << "systematics: " << systName << " , calculator: " << calcName << ", hierarchy: " << hierarchy << std::endl;
66 
67 
68  // quantiles
69  std::string InDirQuant = "/pnfs/nova/persistent/analysis/numu/Ana2018/provisional/quantiles";
70  std::string InFileQuant = InDirQuant + "/quantiles__" + polarity + "_full__numu2018.root";
71  std::cout << "\n\n --- get quantile cuts from file ---" << std::endl;
72  TFile* inFile = TFile::Open( pnfs2xrootd(InFileQuant).c_str() );
73  TH2* FDSpec2D = (TH2*)inFile->FindObjectAny("FDSpec2D");
74  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, 4);
75  inFile->Close();
76  delete inFile;
77 
78 
79  //////////////////// make loader for prediction ////////////////////
80  std::cout << "loaders" << std::endl;
83 
84  std::cout << "set loaders" << std::endl;
85  set_loaders_ana2018(loaders, polarity, period, production);
86 
87 
88  std::cout << "declaring NumuExtrapGenerator and PredictionInterp" << std::endl;
89  std::vector<NumuExtrapGenerator*> NumuCCExtrapVec;
90  std::vector<PredictionInterp*> predictionVec;
91 
92  // hade cuts
93  for(auto thisCut : HadEFracQuantCuts){
94  NumuCCExtrapVec.push_back(new NumuExtrapGenerator(kNumuCCOptimisedAxis, thisCut && kNumuFDCut, thisCut && kNumuNDCut, kNoShift, kPPFXFluxCVWgt*kXSecWeight));
95  predictionVec.push_back(new PredictionInterp(systs, calc, *(NumuCCExtrapVec.back()), loaders));
96  }
97 
98  std::cout << "loaders.Go()" << std::endl;
99  loaders.Go();
100 
101 
102  std::cerr << "Saving Predictions" << std::endl;
103  for(size_t quant=1; quant <= predictionVec.size();quant++){
104  auto prediction = predictionVec[quant-1];
105  std::string OutFileName = Form("prediction_quant%d__", quant) + extrapName + "_" + systName + "__" + polarity + "_" + calcName + "_" + hierarchy + "__" + period + "__numu2018.root";
106  std::cerr << "Saving quantile " << quant << " to file: " << OutFileName << std::endl;
107  TFile *f = TFile::Open(OutFileName.c_str(), "RECREATE");
108  prediction -> SaveTo(f, "prediction");
109  f->Close();
110  }
111 
112 
113 }//make_predictions_extrap
Implements systematic errors by interpolation between shifted templates.
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Cut kNumuCutND2018
Definition: NumuCuts2018.h:41
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
T sqrt(T number)
Definition: d0nt_math.hpp:156
void set_loaders_ana2018(Loaders &loaders, std::string polarity, std::string period, std::string production)
const HistAxis kHadEFracAxis("E_{had.} / E_{#nu}", Binning::Simple(200, 0, 1), kHadEFrac)
HistAxis that implements Hadronic Energy fraction binning used by L.Vinton to derive Hadronic Energy ...
Definition: HistAxes.h:30
OStream cerr
Definition: OStream.cxx:7
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: UtilsExt.cxx:237
osc::OscCalcDumb calc
void make_prediction_extrap(std::string polarity="rhc", std::string period="full", bool nh=true, bool systematics=false, bool calc3a=true)
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
virtual void SetDmsq32(const T &dmsq32)=0
Generates extrapolated Numu predictions.
std::string extrapName
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
ifstream inFile
Definition: AnaPlotMaker.h:34
void SetSpillCut(const SpillCut &cut)
Definition: Loaders.cxx:121
std::vector< const ISyst * > getAllAna2018Systs(const EAnaType2018 ana, const bool smallgenie, const BeamType2018 beam, bool isFit)
std::string OutDir
const HistAxis kNumuCCOptimisedAxis("Reconstructed Neutrino Energy (GeV)", kNumuCCEOptimisedBinning, kCCE)
HistAxis that implements optimised numuCCE from L. Vinton. See docdb 16332. This was close to &#39;custC&#39;...
Definition: HistAxes.h:24
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
const Var kXSecCVWgt2018
Definition: XsecTunes.h:49
std::vector< Loaders * > loaders
Definition: syst_header.h:386
virtual void SetTh23(const T &th23)=0
std::vector< Cut > QuantileCutsFromTH2(TH2 *quantileHist, const HistAxis &independentAxis, const HistAxis &quantileAxis, const unsigned int &numQuantiles, const bool verbose)
: Do the same as the QuantileCuts function but taking in the TH2 instead of making it...
bool systematics
Definition: fnexvscaf.C:31
const Cut kNumuCutFD2018
Definition: NumuCuts2018.h:39
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
T asin(T number)
Definition: d0nt_math.hpp:60
enum BeamMode string