getPredictions.C
Go to the documentation of this file.
3 #include "CAFAna/Core/Loaders.h"
5 #include "CAFAna/Cuts/Cuts.h"
22 #include "CAFAna/Vars/Vars.h"
23 #include "OscLib/OscCalcPMNSOpt.h"
24 
25 using namespace ana;
26 
27 const unsigned int kNumQuantiles = 4;
28 
29 void getPredictions(std::string hornCurrStr, std::string outDir, bool splitNDByQuants = false)
30 {
31  Loaders::FluxType hornCurr;
32  if (hornCurrStr == "FHC")
33  hornCurr = Loaders::kFHC;
34  else if (hornCurrStr == "RHC")
35  hornCurr = Loaders::kRHC;
36  else
37  {
38  std::cerr << "Unrecognized horn current string: '" << hornCurrStr << "'. Options are FHC, RHC." << std::endl;
39  return;
40  }
41  std::transform(hornCurrStr.begin(), hornCurrStr.end(), hornCurrStr.begin(), ::tolower);
42 
43  std::cout << "Making predictions..." << std::endl;
44 
46  // the ND "data" needs to be ND MC for the purposes of the tests in this macro
49 
52 
53  std::vector<const ISyst*> systs = getAllDirectNumuSysts2018();
54 
55  const unsigned int kNumSysts = systs.size();
56 
57  TFile* f = new TFile(Form("%s/Prediction_ccE_%s.root", outDir.c_str(), hornCurrStr.c_str()), "RECREATE");
58 
59  const Cut kNDCut = kNumuCutND2018;
60  const Cut kFDCut = kNumuCutFD2018 && kHasNeutrino;
61  const Var kWeights = kPPFXFluxCVWgt * kXSecCVWgt2018;
62 
63  auto cvmfs_dir = std::getenv("NUMUDATA_DIR");
64  if (!cvmfs_dir)
65  {
66  std::cerr << "Couldn't find UPS dir for numu prediction!" << std::endl;
67  return;
68  }
69  std::string fdspecfile = std::string(cvmfs_dir) + "/ana2018/Quantiles/quantiles__" + hornCurrStr + "_full__numu2018.root";
70 
71  TFile* inFile = TFile::Open( fdspecfile.c_str() );
72  TH2 *FDSpec2D = (TH2*)inFile->FindObjectAny( "FDSpec2D" );
73  std::vector<Cut> kQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, kNumQuantiles);
74  kQuantCuts.emplace_back( kNoCut );
75 
76  PredictionInterp* prediction[kNumSysts*2][kNumQuantiles+1];
77  PredictionNoExtrap* predictionFD[kNumSysts*2][kNumQuantiles+1];
78  PredictionNoExtrap* nominalPred[kNumQuantiles+1];
79 
80  for(unsigned int i = 0; i < kNumQuantiles +1; i++)
81  {
82  nominalPred[i] = new PredictionNoExtrap(loaders, kNumuCCOptimisedAxis, kFDCut && kQuantCuts[i], kNoShift, kWeights);
83 
84  for(unsigned int k = 0; k < kNumSysts; k++)
85  {
86  SystShifts shiftsPos = SystShifts::Nominal();
87  shiftsPos.SetShift(systs[k], 1);
88  NumuExtrapGenerator NumuCCExtrapPos( kNumuCCOptimisedAxis, kFDCut && kQuantCuts[i], splitNDByQuants ? kNDCut && kQuantCuts[i] : kNDCut, shiftsPos, kWeights );
89  prediction[2*k][i] = new PredictionInterp({}, &calc, NumuCCExtrapPos, loaders);
90  predictionFD[2*k][i] = new PredictionNoExtrap(loaders, kNumuCCOptimisedAxis, kFDCut && kQuantCuts[i], shiftsPos, kWeights);
91 
92  SystShifts shiftsNeg = SystShifts::Nominal();
93  shiftsNeg.SetShift(systs[k], -1);
94  NumuExtrapGenerator NumuCCExtrapNeg( kNumuCCOptimisedAxis, kFDCut && kQuantCuts[i], splitNDByQuants ? kNDCut && kQuantCuts[i] : kNDCut, shiftsNeg, kWeights );
95  prediction[2*k+1][i] = new PredictionInterp({}, &calc, NumuCCExtrapNeg, loaders);
96  predictionFD[2*k+1][i] = new PredictionNoExtrap(loaders, kNumuCCOptimisedAxis, kFDCut && kQuantCuts[i], shiftsNeg, kWeights);
97  } // end loop over systs
98  } // end loop over quantiles
99 
100  // GO!
101  loaders.Go();
102 
103  std::string predNameND = "ccE_shiftNDdata";
104  std::string predNameFD = "ccE_noExtrapFD";
105 
106  for(unsigned int i = 0; i < kNumQuantiles +1; i++)
107  {
108  nominalPred[i]->SaveTo( f, Form("nominalPred_%i",i) ) ;
109 
110  for(unsigned int k = 0; k < kNumSysts; k++)
111  {
112  prediction[2*k][i]->SaveTo(f, Form( (predNameND+"_Pos_"+systs[k]->ShortName()+"_%i").c_str(), i) );
113  prediction[2*k+1][i]->SaveTo(f, Form( (predNameND+"_Neg_"+systs[k]->ShortName()+"_%i").c_str(), i) );
114 
115  predictionFD[2*k][i]->SaveTo(f, Form( (predNameFD+"_Pos_"+systs[k]->ShortName()+"_%i").c_str(), i) );
116  predictionFD[2*k+1][i]->SaveTo(f, Form( (predNameFD+"_Neg_"+systs[k]->ShortName()+"_%i").c_str(), i) );
117  } // end loop over systs
118  } // end loop over quantiles
119 
120  f->Close();
121 
122  std::cout << "All predictions made!" << std::endl;
123 
124 } // End of function
Near Detector underground.
Definition: SREnums.h:10
Implements systematic errors by interpolation between shifted templates.
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
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
const Cut kNumuCutND2018
Definition: NumuCuts2018.h:41
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
static SystShifts Nominal()
Definition: SystShifts.h:34
osc::OscCalcDumb calc
std::string outDir
const unsigned int kNumQuantiles
Generates extrapolated Numu predictions.
const Cut kNDCut
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
std::vector< const ISyst * > getAllDirectNumuSysts2018()
Definition: NumuSysts.cxx:114
std::string GetLoaderPath(caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap) const
Definition: Loaders.cxx:89
ifstream inFile
Definition: AnaPlotMaker.h:34
void SetSpillCut(const SpillCut &cut)
Definition: Loaders.cxx:121
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
std::string getenv(std::string const &name)
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
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:22
OStream cout
Definition: OStream.cxx:6
const Cut kHasNeutrino([](const caf::SRProxy *sr){return(sr->mc.nnu!=0);})
Check if MC slice has neutrino information (useful for in-and-out tests)
Definition: TruthCuts.h:61
void getPredictions(std::string hornCurrStr, std::string outDir, bool splitNDByQuants=false)
const Var kXSecCVWgt2018
Definition: XsecTunes.h:48
std::vector< Loaders * > loaders
Definition: syst_header.h:386
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...
const Cut kNumuCutFD2018
Definition: NumuCuts2018.h:39
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
Prediction that just uses FD MC, with no extrapolation.
For nominal spectra and reweighting systs (xsec/flux)
Definition: Prod4Loaders.h:96
void SetLoaderPath(const std::string &path, caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Configure loader via wildcard path.
Definition: Loaders.cxx:25
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:81
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:109
enum BeamMode string