makePred.C
Go to the documentation of this file.
1 #include "CAFAna/Core/Loaders.h"
3 #include "CAFAna/Core/Spectrum.h"
4 
5 #include "CAFAna/Cuts/Cuts.h"
11 
13 #include "CAFAna/Vars/Vars.h"
16 
18 
19 #include "CAFAna/Analysis/Calcs.h"
20 #include "OscLib/OscCalcPMNSOpt.h"
21 
26 
27 #include "TFile.h"
28 
29 
30 using namespace ana;
31 
32 
33 void setLoadersPathNumu2017(Loaders& loaders, bool fakenddata){
34 
35  // std::string nd_data = "prod_sumdecaf_R17-03-01-prod3reco.d_nd_numi_fhc_full_v1_goodruns_numu2017";
36  // std::string nd_nonswap = "prod_sumdecaf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1_numu2017";
37  // std::string fd_nonswap = "prod_sumdecaf_R17-03-01-prod3reco.l_fd_genie_nonswap_fhc_nova_v08_full_v1_numu2017";
38  std::string nd_data = "prod_sumdecaf_R17-03-01-prod3reco.d_nd_numi_fhc_full_v1_goodruns_numu_sa_2017";
39  std::string nd_nonswap = "prod_sumdecaf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1_numu_sa_2017";
40  std::string fd_nonswap = "prod_sumdecaf_R17-03-01-prod3reco.l_fd_genie_nonswap_fhc_nova_v08_full_v1_numu_sa_2017";
41 
42  std::string fd_fluxswap = "prod_sumdecaf_R17-03-01-prod3reco.l_fd_genie_fluxswap_fhc_nova_v08_full_v1_numu2017";
43  std::string fd_tau = "prod_sumdecaf_R17-03-01-prod3reco.l_fd_genie_tau_fhc_nova_v08_full_v1_numu2017";
44 
45  //////////////////// set loaders ////////////////////
46  if(fakenddata) loaders.SetLoaderPath( nd_nonswap, caf::kNEARDET,
49  if(!fakenddata) loaders.SetLoaderPath( nd_data, caf::kNEARDET,
52 
53  loaders.SetLoaderPath( nd_nonswap, caf::kNEARDET, Loaders::kMC, ana::kBeam,
55  loaders.SetLoaderPath( fd_nonswap, caf::kFARDET, Loaders::kMC, ana::kBeam,
57  loaders.SetLoaderPath( fd_fluxswap, caf::kFARDET, Loaders::kMC, ana::kBeam,
61 
62 }
63 
64 void makePred(bool usesysts = true, bool useFakeNDData = false){
65 
66  std::cout << "\n use full numu systs?: " << usesysts
67  << "\n " <<std::endl;
68 
69  // custom osc. parameters (non-max mixing) Get the numu result parameters!
70  osc::OscCalcPMNSOpt calc_nonMax;
71  ana::ResetOscCalcToDefault(&calc_nonMax);
72  calc_nonMax.SetDmsq32(2.67e-3);
73  calc_nonMax.SetTh23(asin(sqrt(0.402))); // PRL-published non-max mixing
74 
75  std::vector<const ISyst*> systs = {};
76  if(usesysts) systs = getAllDirectNumuSysts2017();
77 
78  // Get FD files for HadEFrac vs. neutrino energy distribution
79  std::string quantpath = "/pnfs/nova/persistent/analysis/numu/Ana2017/all_FD_histo_for_quantile_cuts.root";
80  TFile* quantfile = TFile::Open(pnfs2xrootd(quantpath).c_str());
81  TH2* quanthist = (TH2*)quantfile->Get("dir_FDSpec2D/FDSpec2D");
82  std::vector<Cut> QuantCuts = QuantileCutsFromTH2(quanthist, kNumuCCOptimisedAxis, kHadEFracAxis, 4);
83 
86  setLoadersPathNumu2017(loaders, useFakeNDData);
87 
88  //NumuExtrapGenerator NumuCCExtrapAllQuants(kNumuCCOptimisedAxis, kNumuCutFD2017, kNumuCutND2017, kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
90  PredictionInterp predictionAllQuants(systs, &calc_nonMax, NumuCCExtrapAllQuants, loaders);
91 
92  // make a vector with an entry for each quantiles
93  std::vector<NumuExtrapGenerator*> NumuCCExtrapVec;
94  std::vector<PredictionInterp*> predictionVec;
95 
96  for(auto thisCut : QuantCuts){
97  //NumuCCExtrapVec.push_back(new NumuExtrapGenerator(kNumuCCOptimisedAxis, thisCut && kNumuCutFD2017, thisCut && kNumuCutND2017, kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017));
98  NumuCCExtrapVec.push_back(new NumuExtrapGenerator(kNumuCCOptimisedAxis, thisCut && kNumuCutFD2017, thisCut && kNumuCutND2017, kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2018));
99  predictionVec.push_back(new PredictionInterp(systs, &calc_nonMax, *(NumuCCExtrapVec.back()), loaders));
100 
101  } // loop over quantiles
102 
103  ////////////////////
104  // GOOOOOOOO!!!!!
105  loaders.Go();
106  ////////////////////
107 
108  // Now we should save the predictions
109  std::cerr << "Saving Predictions" << std::endl;
110  std::string ndDataType = "";
111  if(useFakeNDData) ndDataType = "fakeNDData_";
112  std::string OutFileName = "pred_3A_prod3_SAor2017decafs_2018Wgt_AllQuants_" + ndDataType + ".root";
113  TFile *fAQ = TFile::Open(OutFileName.c_str(), "RECREATE");
114  predictionAllQuants.SaveTo(fAQ, "prediction");
115  fAQ->Close();
116  delete fAQ;
117 
118  for(size_t quant=1; quant <= predictionVec.size();quant++){
119  auto prediction = predictionVec[quant-1];
120  OutFileName = "pred_3A_prod3_SAor2017decafs_2018Wgt_" + (std::string)Form("Quant%d", quant) + "_" + ndDataType + ".root";
121  std::cerr << "Saving Quantile " << quant << " to Prediction file: " << OutFileName << std::endl;
122  TFile *f = TFile::Open(OutFileName.c_str(), "RECREATE");
123  prediction -> SaveTo(f, "prediction");
124  f->Close();
125  delete f;
126  }
127 
128 }
Near Detector underground.
Definition: SREnums.h:10
Implements systematic errors by interpolation between shifted templates.
Far Detector at Ash River.
Definition: SREnums.h:11
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
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Cut kNumuCutFD2017
Definition: NumuCuts2017.h:39
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
T sqrt(T number)
Definition: d0nt_math.hpp:156
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
Generates extrapolated Numu predictions.
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
void SetSpillCut(const SpillCut &cut)
Definition: Loaders.cxx:121
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
void makePred(bool usesysts=true, bool useFakeNDData=false)
Definition: makePred.C:64
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
void setLoadersPathNumu2017(Loaders &loaders, bool fakenddata)
Definition: makePred.C:33
const Cut kNumuCutND2017
Definition: NumuCuts2017.h:41
void SetTh23(const T &th23) 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:21
OStream cout
Definition: OStream.cxx:6
void SetDmsq32(const T &dmsq32) override
const Var kXSecCVWgt2018
Definition: XsecTunes.h:49
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 SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
std::vector< const ISyst * > getAllDirectNumuSysts2017()
Get a vector of all the numu group systs.
Definition: NumuSysts.cxx:168
Float_t e
Definition: plot.C:35
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
T asin(T number)
Definition: d0nt_math.hpp:60
enum BeamMode string