getPredictions_updatedAna.C
Go to the documentation of this file.
2 #include "CAFAna/Core/Loaders.h"
5 #include "CAFAna/Core/Spectrum.h"
7 #include "CAFAna/Cuts/Cuts.h"
10 #include "CAFAna/Cuts/SpillCuts.h"
11 #include "CAFAna/Cuts/TimingCuts.h"
12 #include "CAFAna/Cuts/TruthCuts.h"
23 #include "CAFAna/Vars/FitVars.h"
28 #include "CAFAna/Vars/Vars.h"
29 #include "OscLib/OscCalcPMNSOpt.h"
30 
31 using namespace ana;
32 
33 const unsigned int kNumQuantiles = 4;
34 
35 void getPredictions_updatedAna(std::string& period, bool splitND = false)
36 {
37 
38  std::cout << "Making predictions..." << std::endl;
39 
40  // std::string nd_data = "prod_sumdecaf_R17-03-01-prod3reco.d_nd_numi_fhc_" + period + "_v1_goodruns_numu2017";
41  std::string nd_nonswap = "prod_sumdecaf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_" + period + "_v1_numu2017";
42  std::string fd_nonswap = "prod_sumdecaf_R17-03-01-prod3reco.l_fd_genie_nonswap_fhc_nova_v08_" + period + "_v1_numu2017";
43  std::string fd_fluxswap = "prod_sumdecaf_R17-03-01-prod3reco.l_fd_genie_fluxswap_fhc_nova_v08_" + period + "_v1_numu2017";
44  std::string fd_tau = "prod_sumdecaf_R17-03-01-prod3reco.l_fd_genie_tau_fhc_nova_v08_" + period + "_v1_numu2017";
45 
47  // loaders.SetLoaderPath( nd_data, caf::kNEARDET, Loaders::kData, ana::kBeam, Loaders::kNonSwap );
54 
57 
58  std::vector<const ISyst*> systs = getAllNumuSystsSA();
59 
60  const unsigned int kNumSysts = systs.size();
61 
62  TFile* f = new TFile(("Predictions/Prediction_" + period + "_ccE.root").c_str(),"RECREATE");
63 
64  const Cut kCVNcut = kCVNm > 0.6;
65  const Cut kBDTcut(
66  [](const caf::SRProxy* sr)
67  { return ( sr->sel.cosrej.numucontpid > 0.5 );
68  }
69  );
70 
71  const Cut kNDCut = kNumuQuality && kCVNcut;
72  const Cut kFDCut = kNumuQuality && kCVNcut && kBDTcut && kHasNeutrino;
73  const Var kWeights = kPPFXFluxCVWgt * kXSecCVWgt2017;
74 
75  SpectrumLoader loader_fd_nonswap(fd_nonswap);
76 
77  std::vector<Cut> kQuantCuts = QuantileCuts(loader_fd_nonswap,
80  kFDCut && kIsNumuCC,
83  kNoShift,
84  kWeights); // currently applying the FD selection to the ND
85 
86  kQuantCuts.push_back(kNoCut);
87 
88  PredictionInterp* prediction[kNumSysts*2][kNumQuantiles+1];
89  PredictionNoExtrap* predictionFD[kNumSysts*2][kNumQuantiles+1];
90  PredictionNoExtrap* nominalPred[kNumQuantiles+1];
91 
92  for(unsigned int i = 0; i < kNumQuantiles +1; i++)
93  {
94  nominalPred[i] = new PredictionNoExtrap(loaders, kNumuCCOptimisedAxis, kFDCut && kQuantCuts[i], kNoShift, kWeights);
95 
96  for(unsigned int k = 0; k < kNumSysts; k++)
97  {
99 
100  shifts.SetShift(systs[k], 1);
101  NumuExtrapGenerator NumuCCExtrapPos( kNumuCCOptimisedAxis, kFDCut && kQuantCuts[i], splitND ? kNDCut && kQuantCuts[i] : kNDCut, shifts, kWeights );
102  prediction[2*k][i] = new PredictionInterp({}, &calc, NumuCCExtrapPos, loaders);
103  predictionFD[2*k][i] = new PredictionNoExtrap(loaders, kNumuCCOptimisedAxis, kFDCut && kQuantCuts[i], shifts, kWeights);
104 
105  shifts.SetShift(systs[k], -1);
106  NumuExtrapGenerator NumuCCExtrapNeg( kNumuCCOptimisedAxis, kFDCut && kQuantCuts[i], splitND ? kNDCut && kQuantCuts[i] : kNDCut, shifts, kWeights );
107  prediction[2*k+1][i] = new PredictionInterp({}, &calc, NumuCCExtrapNeg, loaders);
108  predictionFD[2*k+1][i] = new PredictionNoExtrap(loaders, kNumuCCOptimisedAxis, kFDCut && kQuantCuts[i], shifts, kWeights);
109  } // end loop over systs
110  } // end loop over quantiles
111 
112  // GO!
113  loaders.Go();
114 
115  std::string predNameND = "ccE_shiftNDdata";
116  std::string predNameFD = "ccE_noExtrapFD";
117 
118  for(unsigned int i = 0; i < kNumQuantiles +1; i++)
119  {
120  nominalPred[i]->SaveTo( f, Form("nominalPred_%i",i) ) ;
121 
122  for(unsigned int k = 0; k < kNumSysts; k++)
123  {
124  prediction[2*k][i]->SaveTo(f, Form( (predNameND+"_Pos_"+systs[k]->ShortName()+"_%i").c_str(), i) );
125  prediction[2*k+1][i]->SaveTo(f, Form( (predNameND+"_Neg_"+systs[k]->ShortName()+"_%i").c_str(), i) );
126 
127  predictionFD[2*k][i]->SaveTo(f, Form( (predNameFD+"_Pos_"+systs[k]->ShortName()+"_%i").c_str(), i) );
128  predictionFD[2*k+1][i]->SaveTo(f, Form( (predNameFD+"_Neg_"+systs[k]->ShortName()+"_%i").c_str(), i) );
129  } // end loop over systs
130  } // end loop over quantiles
131 
132  f->Close();
133 
134  std::cout << "All predictions made!" << std::endl;
135 
136 } // End of function
Near Detector underground.
Definition: SREnums.h:10
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
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
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
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
std::vector< Cut > QuantileCuts(SpectrumLoader &loader, const HistAxis &independentAxis, const HistAxis &quantileAxis, const Cut &cut, const unsigned int &numQuantiles, const SpillCut &spillCut, const SystShifts &shift, const Var &wei, const bool verbose)
Returns a vector of Cuts, each one for a different quantile An individual cut will only pass events f...
static SystShifts Nominal()
Definition: SystShifts.h:34
osc::OscCalcDumb calc
Generates extrapolated Numu predictions.
const Cut kNDCut
caf::Proxy< caf::SRCosRej > cosrej
Definition: SRProxy.h:1252
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
const unsigned int kNumQuantiles
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
caf::StandardRecord * sr
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:21
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
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
void getPredictions_updatedAna(std::string &period, bool splitND=false)
std::vector< Loaders * > loaders
Definition: syst_header.h:386
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.
caf::Proxy< caf::SRIDBranch > sel
Definition: SRProxy.h:2141
const Cut kNumuQuality
Definition: NumuCuts.h:18
const Var kXSecCVWgt2017
Definition: XsecTunes.h:37
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
const Var kCVNm
PID
Definition: Vars.cxx:39
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:80
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:109
enum BeamMode string