makePred.C
Go to the documentation of this file.
1 // Makes all the predictions I need:
2 // SA-like, A17-like (per quartile) and OR
3 
4 #include "CAFAna/Core/Loaders.h"
6 #include "CAFAna/Core/Spectrum.h"
7 
8 #include "CAFAna/Cuts/Cuts.h"
12 
14 #include "CAFAna/Vars/Vars.h"
18 
19 #include "CAFAna/Analysis/Calcs.h"
20 #include "OscLib/OscCalcPMNSOpt.h"
21 
25 
26 
27 #include "TFile.h"
28 
29 
30 using namespace ana;
31 
33 
34  std::string nd_data = "prod_sumdecaf_R17-03-01-prod3reco.d_nd_numi_fhc_full_v1_goodruns_numu2017";
35  std::string nd_nonswap = "prod_sumdecaf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1_numu2017";
36  std::string fd_nonswap = "prod_caf_R17-03-01-prod3reco.l_fd_genie_nonswap_fhc_nova_v08_full_v1";
37  std::string fd_fluxswap = "prod_caf_R17-03-01-prod3reco.l_fd_genie_fluxswap_fhc_nova_v08_full_v1";
38  std::string fd_tau = "prod_caf_R17-03-01-prod3reco.l_fd_genie_tau_fhc_nova_v08_full_v1";
39 
40  //////////////////// set loaders ////////////////////
43  loaders.SetLoaderPath( nd_nonswap, caf::kNEARDET, Loaders::kMC, ana::kBeam,
45  loaders.SetLoaderPath( fd_nonswap, caf::kFARDET, Loaders::kMC, ana::kBeam,
47  loaders.SetLoaderPath( fd_fluxswap, caf::kFARDET, Loaders::kMC, ana::kBeam,
51 
52 }
53 
54 void makePred(){
55 
56  // custom osc. parameters (non-max mixing) Get the numu result parameters!
57  osc::OscCalcPMNSOpt calc_nonMax;
58  ana::ResetOscCalcToDefault(&calc_nonMax);
59  calc_nonMax.SetDmsq32(2.43071e-3);
60  calc_nonMax.SetTh23(asin(sqrt(0.471974)));
61  calc_nonMax.SetdCP(M_PI*0.915869); // preliminary Ana2017
62 
63  // Get FD files for HadEFrac vs. neutrino energy distribution
64  std::string quantpath = "/pnfs/nova/persistent/analysis/numu/Ana2017/all_FD_histo_for_quantile_cuts.root";
65  TFile* quantfile = TFile::Open(pnfs2xrootd(quantpath).c_str());
66  TH2* quanthist = (TH2*)quantfile->Get("dir_FDSpec2D/FDSpec2D");
67  std::vector<Cut> QuantCuts = QuantileCutsFromTH2(quanthist, kNumuCCOptimisedAxis, kHadEFracAxis, 4);
68 
71  setLoadersPathNumu2017(loaders);
72 
73  const Cut kNumuCosmicRejSA(
74  [](const caf::SRProxy* sr)
75  { return ( sr->sel.cosrej.anglekal > 0.5 &&
76  sr->sel.cosrej.numuSAcontpid > 0.535 &&
77  //sr->sel.cosrej.numucontpid > 0.535 &&
78  sr->slc.nhit < 400);
79  }
80  );
81 
82  const Cut kNumuFDSA = kNumuQuality && kNumuContainFD && kNumuNCRej && kNumuCosmicRejSA;
83 
84  // make a vector with an entry for each quantiles
85  std::vector<NumuExtrapGenerator*> NumuCCExtrapVec;
86  std::vector<PredictionInterp*> predictionVec;
87 
89  NumuExtrapGenerator NumuCCExtrap_OR (kNumuCCOptimisedAxis, kNumuCutFD2017 || kNumuFDSA, kNumuCutND2017, kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
90 
91  PredictionInterp prediction_SA ({}, &calc_nonMax, NumuCCExtrap_SA, loaders);
92  PredictionInterp prediction_OR ({}, &calc_nonMax, NumuCCExtrap_OR, loaders);
93 
94  for(auto thisCut : QuantCuts){
95  NumuCCExtrapVec.push_back(new NumuExtrapGenerator(kNumuCCOptimisedAxis, thisCut && kNumuCutFD2017, thisCut && kNumuCutND2017, kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017));
96  predictionVec.push_back(new PredictionInterp({}, &calc_nonMax, *(NumuCCExtrapVec.back()), loaders));
97  } // loop over quantiles
98 
99  ////////////////////
100  // GOOOOOOOO!!!!!
101  loaders.Go();
102  ////////////////////
103 
104  // Now we should save the predictions
105  std::cerr << "Saving Predictions" << std::endl;
106  std::string OutFileName = "pred_3A_prod3_AllQuants.root";
107 
108  TFile *f = TFile::Open(OutFileName.c_str(), "RECREATE");
109  prediction_SA.SaveTo(f, "prediction_SA");
110  prediction_OR.SaveTo(f, "prediction_OR");
111  f->Close();
112  delete f;
113 
114 
115  for(size_t quant=1; quant <= predictionVec.size();quant++){
116  auto prediction = predictionVec[quant-1];
117  OutFileName = "pred_3A_prod3_" + (std::string)Form("Quant%d", quant) + ".root";
118  std::cerr << "Saving Quantile " << quant << " to Prediction file: " << OutFileName << std::endl;
119  TFile *f = TFile::Open(OutFileName.c_str(), "RECREATE");
120  prediction -> SaveTo(f, "prediction");
121  f->Close();
122  delete f;
123  }
124 
125 }
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
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
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
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
const Cut kNumuContainFD([](const caf::SRProxy *sr){ std::pair< int, int > planes=calcFirstLastLivePlane(sr->slc.firstplane, std::bitset< 14 >(sr->hdr.dibmask));int planestofront=sr->slc.firstplane-planes.first;int planestoback=planes.second-sr->slc.lastplane;return( sr->slc.ncellsfromedge > 1 &&planestofront > 1 &&planestoback > 1 &&sr->sel.contain.kalfwdcell > 10 &&sr->sel.contain.kalbakcell > 10 &&sr->sel.contain.cosfwdcell > 0 &&sr->sel.contain.cosbakcell > 0);})
Definition: NumuCuts.h:20
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: UtilsExt.cxx:237
#define M_PI
Definition: SbMath.h:34
Generates extrapolated Numu predictions.
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
void makePred(bool usesysts=true, bool useFakeNDData=false)
Definition: makePred.C:64
caf::Proxy< unsigned int > nhit
Definition: SRProxy.h:1315
const HistAxis kNumuCCAxis("Reconstructed Neutrino Energy (GeV)", kNumuEnergyBinning, kCCE)
Definition: HistAxes.h:8
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
caf::StandardRecord * sr
void setLoadersPathNumu2017(Loaders &loaders, bool fakenddata)
Definition: makePred.C:33
const Cut kNumuCutND2017
Definition: NumuCuts2017.h:41
void SetTh23(const T &th23) override
const Cut kNumuNCRej([](const caf::SRProxy *sr){return(sr->sel.remid.pid >0.75);})
Definition: NumuCuts.h:24
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
void SetDmsq32(const T &dmsq32) override
caf::Proxy< caf::SRSlice > slc
Definition: SRProxy.h:2142
void SetdCP(const T &dCP) override
std::vector< Loaders * > loaders
Definition: syst_header.h:386
caf::Proxy< float > anglekal
Definition: SRProxy.h:859
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
caf::Proxy< caf::SRIDBranch > sel
Definition: SRProxy.h:2141
const Cut kNumuQuality
Definition: NumuCuts.h:18
const Var kXSecCVWgt2017
Definition: XsecTunes.h:37
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