make_quantile_cuts_predictions.C
Go to the documentation of this file.
1 #include "CAFAna/Core/Loaders.h"
3 #include "CAFAna/Core/Spectrum.h"
4 
6 #include "CAFAna/Analysis/SALoaders.h"
7 
8 #include "CAFAna/Cuts/Cuts.h"
10 #include "CAFAna/Cuts/TimingCuts.h"
13 
16 
18 #include "CAFAna/Vars/Vars.h"
20 
23 #include "OscLib/OscCalcPMNSOpt.h"
24 
25 #include "TFile.h"
26 
27 using namespace ana;
28 
29 
30 std::vector< std::string > getFDFileVector(std::string kEpoch){
31 
32  std::vector< std::string > fileList;
33  std::string filename = "/pnfs/nova/persistent/production/concat/R16-03-03-prod2reco.f/" + kEpoch + "/prod_decaf_R16-03-03-prod2reco.f_fd_genie_nonswap_fhc_nova_v08_" + kEpoch + "_numu_contain_v1_prod2-snapshot.root";
34 
35  fileList.push_back(filename);
36  return fileList;
37 }
38 
39 
40 void setLoadersPathConcat(SADecafLoaders& loaders, std::string epoch){
41  // feed string from: period{1,2}, epoch{3b,3c} or ""
42 
43  std::string inputDird = "/pnfs/nova/persistent/production/concat/R16-03-03-prod2reco.d/";
44  std::string inputDirf = "/pnfs/nova/persistent/production/concat/R16-03-03-prod2reco.f/";
45 
46  //////////////////// ND, no period splitting ////////////////////
47  std::string nd_data_concat = "prod_decaf_R16-03-03-prod2reco.d_nd_numi_fhc_full_numu_contain_v1_goodruns.root";
48 
49  std::vector< std::string > nd_nonswap_vec;
50  for(int i = 1; i<=16; i++){
51  std::string nd_nonswap_concat = inputDird + "prod_decaf_R16-03-03-prod2reco.d_nd_genie_nonswap_genierw_fhc_nova_v08_full_numu_contain_v1/prod_decaf_R16-03-03-prod2reco.d_nd_genie_nonswap_genierw_fhc_nova_v08_full_numu_contain_v1_" + Form("%d",i) + "_of_16.root";
52  nd_nonswap_vec.push_back(nd_nonswap_concat);
53  }
54 
55  //////////////////// FD, split by period and epoch ////////////////////
56  std::vector< std::string > fd_nonswap_concat_vec;
57  std::vector< std::string > fd_fluxswap_concat_vec;
58  std::vector< std::string > fd_tau_concat_vec;
59 
60  std::string fd_nonswap_concat = epoch + "/prod_decaf_R16-03-03-prod2reco.f_fd_genie_nonswap_fhc_nova_v08_" +
61  epoch + "_numu_contain_v1_prod2-snapshot.root";
62  std::string fd_fluxswap_concat = epoch + "/prod_decaf_R16-03-03-prod2reco.f_fd_genie_fluxswap_fhc_nova_v08_" +
63  epoch + "_numu_contain_v1_prod2-snapshot.root";
64  std::string fd_tau_concat = epoch + "/prod_decaf_R16-03-03-prod2reco.f_fd_genie_tau_fhc_nova_v08_" +
65  epoch + "_numu_contain_v1_prod2-snapshot.root";
66 
67  fd_nonswap_concat_vec.push_back(inputDirf + fd_nonswap_concat);
68  fd_fluxswap_concat_vec.push_back(inputDirf + fd_fluxswap_concat);
69  fd_tau_concat_vec.push_back(inputDirf + fd_tau_concat);
70 
71 
72  //////////////////// set loaders ////////////////////
73  loaders.SetLoaderFiles( {inputDird + nd_data_concat}, caf::kNEARDET, SADecafLoaders::kData, ana::kBeam, SADecafLoaders::kNonSwap );
74  loaders.SetLoaderFiles( nd_nonswap_vec, caf::kNEARDET, SADecafLoaders::kMC, ana::kBeam, SADecafLoaders::kNonSwap );
75  loaders.SetLoaderFiles( fd_nonswap_concat_vec, caf::kFARDET, SADecafLoaders::kMC, ana::kBeam, SADecafLoaders::kNonSwap );
76  loaders.SetLoaderFiles( fd_fluxswap_concat_vec, caf::kFARDET, SADecafLoaders::kMC, ana::kBeam, SADecafLoaders::kFluxSwap);
77  loaders.SetLoaderFiles( fd_tau_concat_vec, caf::kFARDET, SADecafLoaders::kMC, ana::kBeam, SADecafLoaders::kTauSwap );
78 
79 }
80 
81 
83 
84 void make_quantile_cuts_predictions(std::string epochname, bool usesysts=false){
85 
86  //Actually do something
87  std::cout << "make_quantile_cuts" << std::endl;
88 
89  //Get a vector of file names to pass to the spectrum loader
90  auto fdFileVec = getFDFileVector(epochname);
91  for(auto file: fdFileVec) std::cout << "file: " << file << std::endl;
92 
93  SpectrumLoader loader(fdFileVec);
94 
95 
96  std::vector<Cut> HadEFracQuantCuts = QuantileCuts(loader, kNumuCCOptimisedAxis, kHadEFracAxis, kNumuSAKTune && kIsNumuCC, 4, kStandardSpillCuts, kNoShift, kTuftsWeightCC);
97 
98  //////////////////// define an oscillation calculator to use ////////////////////
99  // custom osc. parameters (non-max mixing) Get the numu result parameters!
100  osc::OscCalcPMNSOpt calc_nonMax;
101  calc_nonMax.SetL(810);
102  calc_nonMax.SetRho(0); // No matter effects
103  calc_nonMax.SetDmsq21(7.59e-5);
104  calc_nonMax.SetDmsq32(2.6746e-3);
105  calc_nonMax.SetTh12(.601);
106  calc_nonMax.SetTh13(.1567);
107  calc_nonMax.SetdCP(0);
108  calc_nonMax.SetTh23(0.68696); // non max mixing (ssqth23 = 0.4022)
109 
110  //////////////////// define systematics used ////////////////////
111  std::vector<const ISyst*> systs = {};
112  if(usesysts) systs = getAllDirectNumuSysts2017();
113 
114  //////////////////// make loader for prediction ////////////////////
115  SANumuDecafLoaders loaders_period;
116  loaders_period.SetSpillCut(kStandardSpillCuts);
117  setLoadersPathConcat(loaders_period, epochname);
118 
119  std::vector<NumuExtrapGenerator*> NumuCCExtrapVec;
120  std::vector<PredictionInterp*> predictionVec;
121 
122  for(auto thisCut : HadEFracQuantCuts){
123  NumuCCExtrapVec.push_back(new NumuExtrapGenerator(kNumuCCOptimisedAxis, thisCut && kNumuSAKTune, kNDCut, kNoShift, kTuftsWeightCC));
124  predictionVec.push_back(new PredictionInterp(systs, &calc_nonMax, *(NumuCCExtrapVec.back()), loaders_period));
125  }
126 
127 
128  loaders_period.Go();
129 
130  //Now we should save the predictions
131  std::string SystName = "WithSysts";
132  if(!usesysts) SystName = "StatsOnly";
133 
134  std::cerr << "Saving Predictions" << std::endl;
135  for(size_t quant=1; quant <= predictionVec.size();quant++){
136  auto prediction = predictionVec[quant-1];
137  std::string OutFileName = "Prediction_NonMax_" + SystName + "_HybridSelector_4HadEFracQuantiles_OptimisedCCEBins_" + Form("Quant%d", quant) + "_" + epochname + ".root";
138  std::cerr << "Saving Quantile " << quant << " to Prediction file: " << OutFileName << std::endl;
139  TFile *f = TFile::Open(OutFileName.c_str(), "RECREATE");
140  prediction -> SaveTo(f, "prediction");
141  f->Close();
142  }
143  //Can't be bothered to clear up the mess of new'd objects!
144 }
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 SetTh13(const T &th13) override
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
void SetL(double L) override
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
const Color_t kMC
OStream cerr
Definition: OStream.cxx:7
string filename
Definition: shutoffs.py:106
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...
const Cut kNumuContainND([](const caf::SRProxy *sr){return( sr->trk.kalman.ntracks > sr->trk.kalman.idxremid &&sr->slc.ncellsfromedge > 1 &&sr->slc.firstplane > 1 &&sr->slc.lastplane< 212 &&sr->trk.kalman.tracks[0].start.Z()< 1150 &&( sr->trk.kalman.tracks[0].stop.Z()< 1275 ||sr->sel.contain.kalyposattrans< 55) &&( sr->energy.numu.ndhadcalcatE +sr->energy.numu.ndhadcaltranE)< 0.03 &&sr->sel.contain.kalfwdcellnd > 4 &&sr->sel.contain.kalbakcellnd > 8);})
Definition: NumuCuts.h:22
Generates extrapolated Numu predictions.
const Cut kNDCut
list fileList
Definition: cafExposure.py:30
const Cut kNumuHybridNCRej([](const caf::SRProxy *sr){return(sr->sel.remid.pid > 0.5 &&sr->sel.cvn.numuid > 0.5);})
Definition: NumuCuts.h:28
void make_quantile_cuts_predictions(std::string epochname, bool usesysts=false)
const Cut kNumuSAKTune
Definition: NumuCuts.h:57
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
void SetTh23(const T &th23) override
loader
Definition: demo0.py:10
void SetDmsq21(const T &dmsq21) 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
void SetDmsq32(const T &dmsq32) override
const Var kTuftsWeightCC
Definition: XsecTunes.h:31
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
void SetdCP(const T &dCP) override
std::vector< Loaders * > loaders
Definition: syst_header.h:386
std::vector< std::string > getFDFileVector(std::string kEpoch)
void SetTh12(const T &th12) override
void setLoadersPathConcat(SADecafLoaders &loaders, std::string epoch)
TFile * file
Definition: cellShifts.C:17
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:169
const Cut kNumuQuality
Definition: NumuCuts.h:18
Float_t e
Definition: plot.C:35
void SetRho(double rho) override
enum BeamMode string