make_RHC_WrongSign_NumuQ1.C
Go to the documentation of this file.
1 // ----------------------------------------------------------------
2 // Make the Predictions at ND for numu quantile 1, 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"
8 #include "CAFAna/Core/SpectrumLoader.h"
9 #include "CAFAna/Cuts/Cuts.h"
10 #include "CAFAna/Cuts/TruthCuts.h"
11 #include "CAFAna/Cuts/SpillCuts.h"
12 #include "3FlavorAna/Cuts/NueCuts2018.h"
13 #include "3FlavorAna/Cuts/WrongSignCuts2018.h"
14 #include "3FlavorAna/Cuts/NumuCuts2018.h"
15 #include "3FlavorAna/Cuts/QuantileCuts.h"
16 #include "3FlavorAna/Vars/HistAxes.h"
17 #include "CAFAna/Vars/Vars.h"
18 #include "3FlavorAna/Vars/NumuVars.h"
19 #include "3FlavorAna/Vars/NueVars.h"
20 #include "CAFAna/Vars/GenieWeights.h"
21 #include "CAFAna/Vars/PPFXWeights.h"
22 #include "3FlavorAna/Vars/WrongSignBDTVars.h"
23 
24 #include "CAFAna/Prediction/PredictionNoExtrap.h"
25 #include "CAFAna/Prediction/PredictionNoOsc.h"
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 
40 void make_RHC_WrongSign_NumuQ1()
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);
46  ld_nueRHC.SetSpillCut(kStandardSpillCuts);
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
62  Var wei = kXSecCVWgt2018 * kPPFXFluxCVWgt;
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  // Q1
112  PredictionNoOsc pQ1NumuFSProtonRS(ld_numuRHC,kNumuCCOptimisedAxis,kSelectRSNumuND0 && HadEFracQuantCuts[0],kNoShift,wei);
113  PredictionNoOsc pQ1NumuFSProtonWS(ld_numuRHC,kNumuCCOptimisedAxis,kSelectWSNumuND0 && HadEFracQuantCuts[0],kNoShift,wei);
114  PredictionNoOsc pQ1NumuCVNProngRS(ld_numuRHC,kNumuCCOptimisedAxis,kSelectRSNumuND1 && HadEFracQuantCuts[0],kNoShift,wei);
115  PredictionNoOsc pQ1NumuCVNProngWS(ld_numuRHC,kNumuCCOptimisedAxis,kSelectWSNumuND1 && HadEFracQuantCuts[0],kNoShift,wei);
116  PredictionNoOsc pQ1NumuBDT3RS(ld_numuRHC,kNumuCCOptimisedAxis,kSelectRSNumuND2 && HadEFracQuantCuts[0],kNoShift,wei);
117  PredictionNoOsc pQ1NumuBDT3WS(ld_numuRHC,kNumuCCOptimisedAxis,kSelectWSNumuND2 && HadEFracQuantCuts[0],kNoShift,wei);
118  // Splitting into NC from nu and antinu using spectra for this analysis
119  Spectrum sQ1NumuFSProtonRSNC(ld_numuRHC,kNumuCCOptimisedAxis,
120  kSelectRSNumuND0 && kIsNC && !kIsAntiNu && HadEFracQuantCuts[0],kNoShift,wei);
121  Spectrum sQ1NumuFSProtonRSNCAnti(ld_numuRHC,kNumuCCOptimisedAxis,
122  kSelectRSNumuND0 && kIsNC && kIsAntiNu && HadEFracQuantCuts[0],kNoShift,wei);
123  Spectrum sQ1NumuCVNProngRSNC(ld_numuRHC,kNumuCCOptimisedAxis,
124  kSelectRSNumuND1 && kIsNC && !kIsAntiNu && HadEFracQuantCuts[0],kNoShift,wei);
125  Spectrum sQ1NumuCVNProngRSNCAnti(ld_numuRHC,kNumuCCOptimisedAxis,
126  kSelectRSNumuND1 && kIsNC && kIsAntiNu && HadEFracQuantCuts[0],kNoShift,wei);
127  Spectrum sQ1NumuBDT3RSNC(ld_numuRHC,kNumuCCOptimisedAxis,
128  kSelectRSNumuND2 && kIsNC && !kIsAntiNu && HadEFracQuantCuts[0],kNoShift,wei);
129  Spectrum sQ1NumuBDT3RSNCAnti(ld_numuRHC,kNumuCCOptimisedAxis,
130  kSelectRSNumuND2 && kIsNC && kIsAntiNu && HadEFracQuantCuts[0],kNoShift,wei);
131  Spectrum sQ1NumuFSProtonWSNC(ld_numuRHC,kNumuCCOptimisedAxis,
132  kSelectWSNumuND0 && kIsNC && !kIsAntiNu && HadEFracQuantCuts[0],kNoShift,wei);
133  Spectrum sQ1NumuFSProtonWSNCAnti(ld_numuRHC,kNumuCCOptimisedAxis,
134  kSelectWSNumuND0 && kIsNC && kIsAntiNu && HadEFracQuantCuts[0],kNoShift,wei);
135  Spectrum sQ1NumuCVNProngWSNC(ld_numuRHC,kNumuCCOptimisedAxis,
136  kSelectWSNumuND1 && kIsNC && !kIsAntiNu && HadEFracQuantCuts[0],kNoShift,wei);
137  Spectrum sQ1NumuCVNProngWSNCAnti(ld_numuRHC,kNumuCCOptimisedAxis,
138  kSelectWSNumuND1 && kIsNC && kIsAntiNu && HadEFracQuantCuts[0],kNoShift,wei);
139  Spectrum sQ1NumuBDT3WSNC(ld_numuRHC,kNumuCCOptimisedAxis,
140  kSelectWSNumuND2 && kIsNC && !kIsAntiNu && HadEFracQuantCuts[0],kNoShift,wei);
141  Spectrum sQ1NumuBDT3WSNCAnti(ld_numuRHC,kNumuCCOptimisedAxis,
142  kSelectWSNumuND2 && kIsNC && kIsAntiNu && HadEFracQuantCuts[0],kNoShift,wei);
143 
144  // ------------------------------------------------------------------------------------------------
145  // Numu DATA
146  Spectrum sQ1NumuData_FSProtonRS(ld_dataRHCNumu,kNumuCCOptimisedAxis,kSelectRSNumuND0 && HadEFracQuantCuts[0],kNoShift,kUnweighted);
147  Spectrum sQ1NumuData_FSProtonWS(ld_dataRHCNumu,kNumuCCOptimisedAxis,kSelectWSNumuND0 && HadEFracQuantCuts[0],kNoShift,kUnweighted);
148  Spectrum sQ1NumuData_CVNProngRS(ld_dataRHCNumu,kNumuCCOptimisedAxis,kSelectRSNumuND1 && HadEFracQuantCuts[0],kNoShift,kUnweighted);
149  Spectrum sQ1NumuData_CVNProngWS(ld_dataRHCNumu,kNumuCCOptimisedAxis,kSelectWSNumuND1 && HadEFracQuantCuts[0],kNoShift,kUnweighted);
150  Spectrum sQ1NumuData_BDT3RS(ld_dataRHCNumu,kNumuCCOptimisedAxis,kSelectRSNumuND2 && HadEFracQuantCuts[0],kNoShift,kUnweighted);
151  Spectrum sQ1NumuData_BDT3WS(ld_dataRHCNumu,kNumuCCOptimisedAxis,kSelectWSNumuND2 && HadEFracQuantCuts[0],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("./PredictionsAndDataWrongSignNumuQ1.root","RECREATE");
160  outFile->cd();
161  pQ1NumuFSProtonRS.SaveTo(outFile, "pQ1NumuFSProtonRS");
162  pQ1NumuFSProtonWS.SaveTo(outFile, "pQ1NumuFSProtonWS");
163  pQ1NumuCVNProngRS.SaveTo(outFile, "pQ1NumuCVNProngRS");
164  pQ1NumuCVNProngWS.SaveTo(outFile, "pQ1NumuCVNProngWS");
165  pQ1NumuBDT3RS.SaveTo(outFile, "pQ1NumuBDT3RS");
166  pQ1NumuBDT3WS.SaveTo(outFile, "pQ1NumuBDT3WS");
167  sQ1NumuFSProtonRSNC.SaveTo(outFile, "sQ1NumuFSProtonRSNC");
168  sQ1NumuFSProtonRSNCAnti.SaveTo(outFile, "sQ1NumuFSProtonRSNCAnti");
169  sQ1NumuCVNProngRSNC.SaveTo(outFile, "sQ1NumuCVNProngRSNC");
170  sQ1NumuCVNProngRSNCAnti.SaveTo(outFile, "sQ1NumuCVNProngRSNCAnti");
171  sQ1NumuBDT3RSNC.SaveTo(outFile, "sQ1NumuBDT3RSNC");
172  sQ1NumuBDT3RSNCAnti.SaveTo(outFile, "sQ1NumuBDT3RSNCAnti");
173  sQ1NumuFSProtonWSNC.SaveTo(outFile, "sQ1NumuFSProtonWSNC");
174  sQ1NumuFSProtonWSNCAnti.SaveTo(outFile, "sQ1NumuFSProtonWSNCAnti");
175  sQ1NumuCVNProngWSNC.SaveTo(outFile, "sQ1NumuCVNProngWSNC");
176  sQ1NumuCVNProngWSNCAnti.SaveTo(outFile, "sQ1NumuCVNProngWSNCAnti");
177  sQ1NumuBDT3WSNC.SaveTo(outFile, "sQ1NumuBDT3WSNC");
178  sQ1NumuBDT3WSNCAnti.SaveTo(outFile, "sQ1NumuBDT3WSNCAnti");
179 
180  // DATA
181  sQ1NumuData_FSProtonRS.SaveTo(outFile, "sNumuData_FSProtonRS");
182  sQ1NumuData_FSProtonWS.SaveTo(outFile, "sNumuData_FSProtonWS");
183  sQ1NumuData_CVNProngRS.SaveTo(outFile, "sNumuData_CVNProngRS");
184  sQ1NumuData_CVNProngWS.SaveTo(outFile, "sNumuData_CVNProngWS");
185  sQ1NumuData_BDT3RS.SaveTo(outFile, "sNumuData_BDT3RS");
186  sQ1NumuData_BDT3WS.SaveTo(outFile, "sNumuData_BDT3WS");
187 
188  outFile->Close();
189 }
quantiles
std::vector< double > quantiles(TH1D *h)
Definition: absCal.cxx:528
ana
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
ana::kPPFXFluxCVWgt
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
ana::kNumuCutND2018
const Cut kNumuCutND2018
Definition: NumuCuts2018.h:41
caf::Proxy< caf::StandardRecord >
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
ana::kAntiNueBDTCVN
const Var kAntiNueBDTCVN
Definition: WrongSignBDTVars.h:43
make_RHC_WrongSign_NumuQ1
void make_RHC_WrongSign_NumuQ1()
Definition: make_RHC_WrongSign_NumuQ1.C:40
ana::SpectrumLoaderBase::SetSpillCut
void SetSpillCut(const SpillCut &cut)
Definition: SpectrumLoaderBase.cxx:64
check_time_usage.float
float
Definition: check_time_usage.py:6
ana::kNueEnergy2018
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
ana::Energy
const Cut Energy
Definition: NCPi0SingleProngCuts.h:313
sr
caf::StandardRecord * sr
Definition: FillTrainingTrees.h:18
ana::kCVNFSProtonScore2018
const Var kCVNFSProtonScore2018([](const caf::SRProxy *sr){return CVNFinalStateScore2018(sr, 2212);})
Proton score from CVN Final State labels, using Prod3Train CVN.
Definition: CVNFinalStates.h:22
nova
Hold drift constants from current run.
Definition: DriftCache.h:17
ana::kCVNSSe
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
ana::kCVNProngProtonCtNumu
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;})
Definition: WrongSignCuts2018.h:17
ana::SpectrumLoader
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
Definition: SpectrumLoader.h:21
ana::kAntiNumuBDTCVN
const Var kAntiNumuBDTCVN
Definition: WrongSignBDTVars.h:55
ana::kXSecCVWgt2018
const Var kXSecCVWgt2018
Definition: XsecTunes.h:48
ana::assert
assert(nhit_max >=nhit_nbins)
ana::kStandardSpillCuts
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
ana::_Cut< caf::SRProxy >
ana::_Var< caf::SRProxy >
string
enum BeamMode string
Generated on Sat Jan 16 2021 14:12:18 for NOvA by   doxygen 1.8.11