make_RHC_WrongSign_Numu_MC.C
Go to the documentation of this file.
1 // ----------------------------------------------------------------
2 // Make the Predictions at ND for numu, MC only,
3 // to be combined with Matt's neutron study
4 // See docDB-29919 for tech note on WS study
5 // ----------------------------------------------------------------
6 
7 #include "CAFAna/Core/Binning.h"
8 #include "CAFAna/Core/Spectrum.h"
10 #include "CAFAna/Cuts/Cuts.h"
11 #include "CAFAna/Cuts/TruthCuts.h"
12 #include "CAFAna/Cuts/SpillCuts.h"
17 #include "CAFAna/Vars/Vars.h"
23 
26 #include "OscLib/OscCalcPMNSOpt.h"
27 
28 #include "TCanvas.h"
29 #include "TH2.h"
30 #include "TLegend.h"
31 #include "TLatex.h"
32 #include "TFile.h"
33 #include "TProfile.h"
34 #include "TString.h"
35 
36 using namespace ana;
37 
39 {
40  // ND MC data set, Prod 4 RHC
41  // nue
42  std::string fname_numuRHC = "prod_sumdecaf_R17-11-14-prod4reco.e_nd_genie_nonswap_rhc_nova_v08_full_v1_numu2018";
43  SpectrumLoader ld_numuRHC(fname_numuRHC);
44  ld_numuRHC.SetSpillCut(kStandardSpillCuts);
45 
46  // Standard weights
48 
49  // Extra Cuts and Vars for analysis
50  // -------------------------------------------------------------------------------------------------
51  const Cut kCutHighCVNSSe = kCVNSSe >= 0.98; // Ana 2018 Hi CVN bin (docDB-27317)
52  const Cut kCutLowCVNSSe = kCVNSSe >= 0.89; // Ana 2018 Lo CVN bin has this as low edge
53 
54  const Cut kSelectRSNumuND0 = kNumuCutND2018 && kCVNFSProtonScore2018<0.49;
55  const Cut kSelectRSNumuND1 = kNumuCutND2018 && !kCVNProngProtonCtNumu;
56  const Cut kSelectRSNumuND2 = kNumuCutND2018 && kAntiNumuBDTCVN>0.11;
57 
58  const Cut kSelectWSNumuND0 = kNumuCutND2018 && kCVNFSProtonScore2018>=0.49;
59  const Cut kSelectWSNumuND1 = kNumuCutND2018 && kCVNProngProtonCtNumu;
60  const Cut kSelectWSNumuND2 = kNumuCutND2018 && kAntiNumuBDTCVN<=0.11;
61 
62  // Nue-style Binnings for BDT
63  const Var kNueBDT3Bin([](const caf::SRProxy *sr)
64  {
65  if( kAntiNueBDTCVN(sr)>0. ) return float(0.5);
66  else return float(1.5);
67 
68  return float(-1.5);
69  });
70 
71  const Var kNueBDT3AnaBinning([&](const caf::SRProxy *sr)
72  {
73  int selBin = kNueBDT3Bin(sr);
74 
75  float nuE = kNueEnergy2018(sr);
76  int nuEBin = nuE/0.5;
77  assert(nuEBin <= 8 && "An event with nuE > 4.5 should never happen");
78 
79  int anaBin = 9*selBin + nuEBin;
80  return anaBin;
81  });
82 
83  const HistAxis kAxisFSProton("NuE Energy / Analysis Bin",Binning::Simple(18,0,18),kNueFSProtonAnaBinning);
84  const HistAxis kAxisCVNProng("NuE Energy / Analysis Bin",Binning::Simple(18,0,18),kNueCVNProngAnaBinning);
85  const HistAxis kAxisBDT3("NuE Energy / Analysis Bin",Binning::Simple(18,0,18),kNueBDT3AnaBinning);
86  // -------------------------------------------------------------------------------------------------
87 
88  // Predictions/Spectra
89  PredictionNoOsc pNumuFSProtonRS(ld_numuRHC,kNumuCCOptimisedAxis,kSelectRSNumuND0,kNoShift,wei);
90  PredictionNoOsc pNumuFSProtonWS(ld_numuRHC,kNumuCCOptimisedAxis,kSelectWSNumuND0,kNoShift,wei);
91  PredictionNoOsc pNumuTot(ld_numuRHC,kNumuCCOptimisedAxis,kNumuCutND2018,kNoShift,wei);
92 
93  //--- Splitting into NC from nu and antinu using spectra for this analysis
94  Spectrum sNumuFSProtonRSNC(ld_numuRHC,kNumuCCOptimisedAxis,
95  kSelectRSNumuND0 && kIsNC && !kIsAntiNu,kNoShift,wei);
96  Spectrum sNumuFSProtonRSNCAnti(ld_numuRHC,kNumuCCOptimisedAxis,
97  kSelectRSNumuND0 && kIsNC && kIsAntiNu,kNoShift,wei);
98  Spectrum sNumuFSProtonWSNC(ld_numuRHC,kNumuCCOptimisedAxis,
99  kSelectWSNumuND0 && kIsNC && !kIsAntiNu,kNoShift,wei);
100  Spectrum sNumuFSProtonWSNCAnti(ld_numuRHC,kNumuCCOptimisedAxis,
101  kSelectWSNumuND0 && kIsNC && kIsAntiNu,kNoShift,wei);
102 
103  Spectrum sNumuTotNC(ld_numuRHC,kNumuCCOptimisedAxis,
104  kNumuCutND2018 && kIsNC && !kIsAntiNu,kNoShift,wei);
105  Spectrum sNumuTotNCAnti(ld_numuRHC,kNumuCCOptimisedAxis,
106  kNumuCutND2018 && kIsNC && kIsAntiNu,kNoShift,wei);
107 
108  // Press play
109  ld_numuRHC.Go();
110 
111  // Save the predictions
112  TFile *outFile = new TFile("./WrongSign-XCheck2019-NumuMC.root","RECREATE");
113  outFile->cd();
114 
115  pNumuFSProtonRS.SaveTo(outFile, "pNumuFSProtonRS");
116  pNumuFSProtonWS.SaveTo(outFile, "pNumuFSProtonWS");
117  pNumuTot.SaveTo(outFile, "pNumuTot");
118 
119  sNumuFSProtonRSNC.SaveTo(outFile, "sNumuFSProtonRSNC");
120  sNumuFSProtonRSNCAnti.SaveTo(outFile, "sNumuFSProtonRSNCAnti");
121  sNumuFSProtonWSNC.SaveTo(outFile, "sNumuFSProtonWSNC");
122  sNumuFSProtonWSNCAnti.SaveTo(outFile, "sNumuFSProtonWSNCAnti");
123  sNumuTotNC.SaveTo(outFile, "sNumuTotNC");
124  sNumuTotNCAnti.SaveTo(outFile, "sNumuTotNCAnti");
125 
126  outFile->Close();
127 }
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Cut kNumuCutND2018
Definition: NumuCuts2018.h:41
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2108
const Var kAntiNueBDTCVN
void SetSpillCut(const SpillCut &cut)
const Var kNueEnergy2018([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;if(kIsRHC(sr)) return kNueEnergyRHC(sr);else return kNueEnergyFHC(sr);})
Definition: NueEnergy2018.h:25
const Cut Energy
caf::StandardRecord * sr
const Var kCVNFSProtonScore2018([](const caf::SRProxy *sr){return CVNFinalStateScore2018(sr, 2212);})
Proton score from CVN Final State labels, using Prod3Train CVN.
const Var kCVNSSe([](const caf::SRProxy *sr){throw std::runtime_error("kCVNSSe is no longer available. Fix your macro so you don't use it.");return-5.;})
2018 nue PID
Definition: Vars.h:52
const Cut kCVNProngProtonCtNumu([](const caf::SRProxy *sr){ int count=0;if(sr->vtx.elastic.fuzzyk.npng==1) return false;for(unsigned int i=0;i< sr->vtx.elastic.fuzzyk.npng;++i){if(sr->vtx.elastic.fuzzyk.png[i].cvnpart.pdgmax!=2212) continue;if(util::pythag(sr->vtx.elastic.vtx.X()-sr->vtx.elastic.fuzzyk.png[i].start.X(), sr->vtx.elastic.vtx.Y()-sr->vtx.elastic.fuzzyk.png[i].start.Y(), sr->vtx.elastic.vtx.Z()-sr->vtx.elastic.fuzzyk.png[i].start.Z()) > 20) continue;if(sr->vtx.elastic.fuzzyk.png[i].cvnpart.maxval< 0.72) continue;count++;}return count >0;})
void make_RHC_WrongSign_Numu_MC()
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Var kAntiNumuBDTCVN
const Var kXSecCVWgt2018
Definition: XsecTunes.h:49
assert(nhit_max >=nhit_nbins)
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49