make_RHC_WrongSign_NumuQ2.C
Go to the documentation of this file.
1 // ----------------------------------------------------------------
2 // Make the Predictions at ND for numu quantile 2, including data
3 // See docDB-29919 for tech note on WS study
4 // ----------------------------------------------------------------
5 
6 #include "CAFAna/Core/Binning.h"
7 #include "CAFAna/Core/Spectrum.h"
9 #include "CAFAna/Cuts/Cuts.h"
10 #include "CAFAna/Cuts/TruthCuts.h"
11 #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 #include "Utilities/rootlogon.C"
37 
38 using namespace ana;
39 
41 {
42  // ND MC data set, Prod 4 RHC
43  // nue
44  std::string fname_nueRHC = "prod_sumdecaf_R17-11-14-prod4reco.e_nd_genie_nonswap_rhc_nova_v08_full_v1_nue2018";
45  SpectrumLoader ld_nueRHC(fname_nueRHC);
47  // numu
48  std::string fname_numuRHC = "prod_sumdecaf_R17-11-14-prod4reco.e_nd_genie_nonswap_rhc_nova_v08_full_v1_numu2018";
49  SpectrumLoader ld_numuRHC(fname_numuRHC);
50  ld_numuRHC.SetSpillCut(kStandardSpillCuts);
51 
52  // ND real data set, all periods concats (I think this is periods 4, 6)
53  std::string rhcDataNue = "prod_sumdecaf_R17-09-05-prod4recopreview.f_nd_numi_rhc_full_v1_addShortSimpleCVN_goodruns_nue2018";
54  SpectrumLoader ld_dataRHCNue(rhcDataNue);
55  ld_dataRHCNue.SetSpillCut(kStandardSpillCuts);
56 
57  std::string rhcDataNumu = "prod_sumdecaf_R17-09-05-prod4recopreview.f_nd_numi_rhc_full_v1_addShortSimpleCVN_goodruns_numu2018";
58  SpectrumLoader ld_dataRHCNumu(rhcDataNumu);
59  ld_dataRHCNumu.SetSpillCut(kStandardSpillCuts);
60 
61  // Standard weights
63 
64  // Extra Cuts and Vars for analysis
65  // -------------------------------------------------------------------------------------------------
66  const Cut kCutHighCVNSSe = kCVNSSe >= 0.98; // Ana 2018 Hi CVN bin (docDB-27317)
67  const Cut kCutLowCVNSSe = kCVNSSe >= 0.89; // Ana 2018 Lo CVN bin has this as low edge
68 
69  const Cut kSelectRSNumuND0 = kNumuCutND2018 && kCVNFSProtonScore2018<0.49;
70  const Cut kSelectRSNumuND1 = kNumuCutND2018 && !kCVNProngProtonCtNumu;
71  const Cut kSelectRSNumuND2 = kNumuCutND2018 && kAntiNumuBDTCVN>0.11;
72 
73  const Cut kSelectWSNumuND0 = kNumuCutND2018 && kCVNFSProtonScore2018>=0.49;
74  const Cut kSelectWSNumuND1 = kNumuCutND2018 && kCVNProngProtonCtNumu;
75  const Cut kSelectWSNumuND2 = kNumuCutND2018 && kAntiNumuBDTCVN<=0.11;
76 
77  // Nue-style Binnings for BDT
78  const Var kNueBDT3Bin([](const caf::SRProxy *sr)
79  {
80  if( kAntiNueBDTCVN(sr)>0. ) return float(0.5);
81  else return float(1.5);
82 
83  return float(-1.5);
84  });
85 
86  const Var kNueBDT3AnaBinning([&](const caf::SRProxy *sr)
87  {
88  int selBin = kNueBDT3Bin(sr);
89 
90  float nuE = kNueEnergy2018(sr);
91  int nuEBin = nuE/0.5;
92  assert(nuEBin <= 8 && "An event with nuE > 4.5 should never happen");
93 
94  int anaBin = 9*selBin + nuEBin;
95  return anaBin;
96  });
97 
98  const HistAxis kAxisFSProton("NuE Energy / Analysis Bin",Binning::Simple(18,0,18),kNueFSProtonAnaBinning);
99  const HistAxis kAxisCVNProng("NuE Energy / Analysis Bin",Binning::Simple(18,0,18),kNueCVNProngAnaBinning);
100  const HistAxis kAxisBDT3("NuE Energy / Analysis Bin",Binning::Simple(18,0,18),kNueBDT3AnaBinning);
101  // -------------------------------------------------------------------------------------------------
102 
103  // Predictions for numu RHC
104  // Numu quantiles
105  TFile *numuCutFile =
106  TFile::Open(pnfs2xrootd("/pnfs/nova/persistent/analysis/numu/Ana2018/provisional/quantiles/quantiles__rhc_full__numu2018.root").c_str());
107  const int NHadEFracQuantiles = 4; // defines how many divisions of hadEFrac are used
108  TH2 *FDSpec2D = (TH2*)numuCutFile->FindObjectAny("FDSpec2D");
109  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, NHadEFracQuantiles);
110 
111  // Q2
112  PredictionNoOsc pQ2NumuFSProtonRS(ld_numuRHC,kNumuCCOptimisedAxis,kSelectRSNumuND0 && HadEFracQuantCuts[1],kNoShift,wei);
113  PredictionNoOsc pQ2NumuFSProtonWS(ld_numuRHC,kNumuCCOptimisedAxis,kSelectWSNumuND0 && HadEFracQuantCuts[1],kNoShift,wei);
114  PredictionNoOsc pQ2NumuCVNProngRS(ld_numuRHC,kNumuCCOptimisedAxis,kSelectRSNumuND1 && HadEFracQuantCuts[1],kNoShift,wei);
115  PredictionNoOsc pQ2NumuCVNProngWS(ld_numuRHC,kNumuCCOptimisedAxis,kSelectWSNumuND1 && HadEFracQuantCuts[1],kNoShift,wei);
116  PredictionNoOsc pQ2NumuBDT3RS(ld_numuRHC,kNumuCCOptimisedAxis,kSelectRSNumuND2 && HadEFracQuantCuts[1],kNoShift,wei);
117  PredictionNoOsc pQ2NumuBDT3WS(ld_numuRHC,kNumuCCOptimisedAxis,kSelectWSNumuND2 && HadEFracQuantCuts[1],kNoShift,wei);
118  // Splitting into NC from nu and antinu using spectra for this analysis
119  Spectrum sQ2NumuFSProtonRSNC(ld_numuRHC,kNumuCCOptimisedAxis,
120  kSelectRSNumuND0 && kIsNC && !kIsAntiNu && HadEFracQuantCuts[1],kNoShift,wei);
121  Spectrum sQ2NumuFSProtonRSNCAnti(ld_numuRHC,kNumuCCOptimisedAxis,
122  kSelectRSNumuND0 && kIsNC && kIsAntiNu && HadEFracQuantCuts[1],kNoShift,wei);
123  Spectrum sQ2NumuCVNProngRSNC(ld_numuRHC,kNumuCCOptimisedAxis,
124  kSelectRSNumuND1 && kIsNC && !kIsAntiNu && HadEFracQuantCuts[1],kNoShift,wei);
125  Spectrum sQ2NumuCVNProngRSNCAnti(ld_numuRHC,kNumuCCOptimisedAxis,
126  kSelectRSNumuND1 && kIsNC && kIsAntiNu && HadEFracQuantCuts[1],kNoShift,wei);
127  Spectrum sQ2NumuBDT3RSNC(ld_numuRHC,kNumuCCOptimisedAxis,
128  kSelectRSNumuND2 && kIsNC && !kIsAntiNu && HadEFracQuantCuts[1],kNoShift,wei);
129  Spectrum sQ2NumuBDT3RSNCAnti(ld_numuRHC,kNumuCCOptimisedAxis,
130  kSelectRSNumuND2 && kIsNC && kIsAntiNu && HadEFracQuantCuts[1],kNoShift,wei);
131  Spectrum sQ2NumuFSProtonWSNC(ld_numuRHC,kNumuCCOptimisedAxis,
132  kSelectWSNumuND0 && kIsNC && !kIsAntiNu && HadEFracQuantCuts[1],kNoShift,wei);
133  Spectrum sQ2NumuFSProtonWSNCAnti(ld_numuRHC,kNumuCCOptimisedAxis,
134  kSelectWSNumuND0 && kIsNC && kIsAntiNu && HadEFracQuantCuts[1],kNoShift,wei);
135  Spectrum sQ2NumuCVNProngWSNC(ld_numuRHC,kNumuCCOptimisedAxis,
136  kSelectWSNumuND1 && kIsNC && !kIsAntiNu && HadEFracQuantCuts[1],kNoShift,wei);
137  Spectrum sQ2NumuCVNProngWSNCAnti(ld_numuRHC,kNumuCCOptimisedAxis,
138  kSelectWSNumuND1 && kIsNC && kIsAntiNu && HadEFracQuantCuts[1],kNoShift,wei);
139  Spectrum sQ2NumuBDT3WSNC(ld_numuRHC,kNumuCCOptimisedAxis,
140  kSelectWSNumuND2 && kIsNC && !kIsAntiNu && HadEFracQuantCuts[1],kNoShift,wei);
141  Spectrum sQ2NumuBDT3WSNCAnti(ld_numuRHC,kNumuCCOptimisedAxis,
142  kSelectWSNumuND2 && kIsNC && kIsAntiNu && HadEFracQuantCuts[1],kNoShift,wei);
143 
144  // ------------------------------------------------------------------------------------------------
145  // Numu DATA
146  Spectrum sQ2NumuData_FSProtonRS(ld_dataRHCNumu,kNumuCCOptimisedAxis,kSelectRSNumuND0 && HadEFracQuantCuts[1],kNoShift,kUnweighted);
147  Spectrum sQ2NumuData_FSProtonWS(ld_dataRHCNumu,kNumuCCOptimisedAxis,kSelectWSNumuND0 && HadEFracQuantCuts[1],kNoShift,kUnweighted);
148  Spectrum sQ2NumuData_CVNProngRS(ld_dataRHCNumu,kNumuCCOptimisedAxis,kSelectRSNumuND1 && HadEFracQuantCuts[1],kNoShift,kUnweighted);
149  Spectrum sQ2NumuData_CVNProngWS(ld_dataRHCNumu,kNumuCCOptimisedAxis,kSelectWSNumuND1 && HadEFracQuantCuts[1],kNoShift,kUnweighted);
150  Spectrum sQ2NumuData_BDT3RS(ld_dataRHCNumu,kNumuCCOptimisedAxis,kSelectRSNumuND2 && HadEFracQuantCuts[1],kNoShift,kUnweighted);
151  Spectrum sQ2NumuData_BDT3WS(ld_dataRHCNumu,kNumuCCOptimisedAxis,kSelectWSNumuND2 && HadEFracQuantCuts[1],kNoShift,kUnweighted);
152  // ------------------------------------------------------------------------------------------------
153 
154  // Press play
155  ld_numuRHC.Go();
156  ld_dataRHCNumu.Go();
157 
158  // Save the predictions
159  TFile *outFile = new TFile("./PredictionsAndDataWrongSignNumuQ2.root","RECREATE");
160  outFile->cd();
161  pQ2NumuFSProtonRS.SaveTo(outFile, "pQ2NumuFSProtonRS");
162  pQ2NumuFSProtonWS.SaveTo(outFile, "pQ2NumuFSProtonWS");
163  pQ2NumuCVNProngRS.SaveTo(outFile, "pQ2NumuCVNProngRS");
164  pQ2NumuCVNProngWS.SaveTo(outFile, "pQ2NumuCVNProngWS");
165  pQ2NumuBDT3RS.SaveTo(outFile, "pQ2NumuBDT3RS");
166  pQ2NumuBDT3WS.SaveTo(outFile, "pQ2NumuBDT3WS");
167  sQ2NumuFSProtonRSNC.SaveTo(outFile, "sQ2NumuFSProtonRSNC");
168  sQ2NumuFSProtonRSNCAnti.SaveTo(outFile, "sQ2NumuFSProtonRSNCAnti");
169  sQ2NumuCVNProngRSNC.SaveTo(outFile, "sQ2NumuCVNProngRSNC");
170  sQ2NumuCVNProngRSNCAnti.SaveTo(outFile, "sQ2NumuCVNProngRSNCAnti");
171  sQ2NumuBDT3RSNC.SaveTo(outFile, "sQ2NumuBDT3RSNC");
172  sQ2NumuBDT3RSNCAnti.SaveTo(outFile, "sQ2NumuBDT3RSNCAnti");
173  sQ2NumuFSProtonWSNC.SaveTo(outFile, "sQ2NumuFSProtonWSNC");
174  sQ2NumuFSProtonWSNCAnti.SaveTo(outFile, "sQ2NumuFSProtonWSNCAnti");
175  sQ2NumuCVNProngWSNC.SaveTo(outFile, "sQ2NumuCVNProngWSNC");
176  sQ2NumuCVNProngWSNCAnti.SaveTo(outFile, "sQ2NumuCVNProngWSNCAnti");
177  sQ2NumuBDT3WSNC.SaveTo(outFile, "sQ2NumuBDT3WSNC");
178  sQ2NumuBDT3WSNCAnti.SaveTo(outFile, "sQ2NumuBDT3WSNCAnti");
179 
180  // DATA
181  sQ2NumuData_FSProtonRS.SaveTo(outFile, "sNumuData_FSProtonRS");
182  sQ2NumuData_FSProtonWS.SaveTo(outFile, "sNumuData_FSProtonWS");
183  sQ2NumuData_CVNProngRS.SaveTo(outFile, "sNumuData_CVNProngRS");
184  sQ2NumuData_CVNProngWS.SaveTo(outFile, "sNumuData_CVNProngWS");
185  sQ2NumuData_BDT3RS.SaveTo(outFile, "sNumuData_BDT3RS");
186  sQ2NumuData_BDT3WS.SaveTo(outFile, "sNumuData_BDT3WS");
187 
188  outFile->Close();
189 }
std::vector< double > quantiles(TH1D *h)
Definition: absCal.cxx:528
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.
Hold drift constants from current run.
Definition: DriftCache.h:17
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_NumuQ2()
::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