test_micheldecomp.C
Go to the documentation of this file.
1 #include "OscLib/OscCalcDumb.h"
2 
3 #include "CAFAna/Vars/Vars.h"
6 
7 #include "CAFAna/Cuts/Cuts.h"
10 #include "CAFAna/Cuts/TruthCuts.h"
11 #include "CAFAna/Cuts/NueCutsSecondAna.h"
12 #include "CAFAna/Core/Spectrum.h"
14 
15 #include "CAFAna/Analysis/SALoaders.h"
25 
26 #include "CAFAna/Systs/Systs.h"
27 #include "CAFAna/Core/SystShifts.h"
28 
29 
30 #include "Utilities/func/MathUtil.h"
31 
32 #include "TCanvas.h"
33 #include "TH2.h"
34 #include "TFile.h"
35 #include "TLegend.h"
36 #include "TProfile.h"
37 #include "TRandom3.h"
39 
40 #include <iostream>
41 #include <cmath>
42 
43 using namespace ana;
44 
46 {
47  if (PID == "LID")
48  return true;
49  if (PID == "LEM")
50  return true;
51  if (PID == "CVN")
52  return true;
53  return false;
54 }
55 
58  return pred.Predict(&calc) - pred.PredictComponent(&calc,
61  Sign::kBoth);
62 }
63 
64 double GetRes(TH1* h1, TH1* h2, int j)
65 {
66  double cont1 = h1->GetBinContent(j);
67  double cont2 = h2->GetBinContent(j);
68 
69  return (cont2 == 0) ? 0 : cont1/cont2 - 1;
70 }
71 
73 {
74  if (PID == "CVN") return kNueAxisSecondAnaCVN;
75  if (PID == "LID") return kNueAxisSecondAnaLID;
76  return kNueAxisSecondAnaLEM;
77 }
78 
80 {
81  if (PID == "CVN") return kNueSecondAnaCVNeSsb;
82  if (PID == "LID") return kNueSecondAnaLIDSsb;
83  return kNueSecondAnaLEMSsb;
84 }
85 
86 void SaveDCMPPlots(const IDecomp* dcmp, std::string label, double pot)
87 {
88  (dcmp->NumuComponent()+dcmp->AntiNumuComponent())
89  .ToTH1(pot)->Write((label+"_cc").c_str());
90  (dcmp->NueComponent()+dcmp->AntiNueComponent())
91  .ToTH1(pot)->Write((label+"_nue").c_str());
92  dcmp->NCTotalComponent().ToTH1(pot)->Write((label+"_nc").c_str());
93 }
94 
96 {
97 
98  assert(CheckPIDString(PID) && "PID must be CVN, LID, or LEM");
99 
100  HistAxis DCMPAxis = GetDCMPAxis(PID);
101  HistAxis NumuAxis("CCE", Binning::Simple(50,0,5), kCCE);
102 
103  Cut PIDCut = GetPIDCut(PID);
104  Cut NDNue = kNueNDSecondAnaPresel && PIDCut;
105  Cut FDNue = kNueSecondAnaFullPresel && PIDCut;
106  Var MCWei = kTuftsWeightCC;
107 
108  SANueExtrapLoaders loaders;
109  loaders.SetSpillCut(kStandardSpillCuts);
110 
111 
112  // Need a numu decomp for all extraps
113  NumuDecomp numu(loaders, NumuAxis, NDNue, kNoShift, kNoShift, MCWei);
114  // Set up MC
115  const CheatDecomp* cheat =
116  new CheatDecomp(loaders.GetLoader(caf::kNEARDET, Loaders::kMC),
117  DCMPAxis, NDNue, kNoShift, MCWei);
118  NueExtrap mc_ex (loaders, *cheat, numu, DCMPAxis,
119  NumuAxis, FDNue, NDNue,kNumuND, kNoShift, MCWei);
120  PredictionExtrap mc_pr(&mc_ex);
121 
122  // Set up MichelDecomp, pegs nues to mc estimate
123  MichelDecomp michel(loaders, DCMPAxis, NDNue, cheat,
124  kNoShift, kNoShift, MCWei);
125  NueExtrap(loaders, michel, numu, DCMPAxis, NumuAxis,
126  FDNue, NDNue,kNumuND, kNoShift, MCWei);
127  PredictionExtrap michel_pr(&michel_ex);
128 
129  // Set up (Michel+BEN)Decomp, pegs nues to estimate from BEN Decomp
130  const BENDecomp* ben = new BENDecomp(loaders, DCMPAxis, NDNue,
131  kNoShift, kNoShift, MCWei);
132  MichelDecomp michelben(loaders, DCMPAxis, NDNue, ben,
133  kNoShift, kNoShift, MCWei);
134  NueExtrap michelben_ex(loaders, michelben, numu, DCMPAxis, NumuAxis,
135  FDNue, NDNue,kNumuND, kNoShift, MCWei);
136  PredictionExtrap michelben_pr(&michelben_ex);
137 
138  // Set up Prop
139  ProportionalDecomp prop(loaders, DCMPAxis, NDNue,
140  kNoShift, kNoShift, MCWei);
141  NueExtrap prop_ex(loaders, prop, numu, DCMPAxis, NumuAxis,
142  FDNue, NDNue, kNumuND, kNoShift, MCWei);
143  PredictionExtrap prop_pr(&prop_ex);
144 
145 
146  loaders.Go();
147 
148 
149  double ndpot = michel.NumuComponent().POT();
150  double fdpot = 6.754e20; // From Evan's DocDB 15089
151 
152  TFile* outFile = new TFile(("mdcmp_"+PID+".root").c_str(),"RECREATE");
153  outFile->cd();
154 
155  BackPredict(mc_pr). ToTH1(fdpot)->Write("FD_mcestimate");
156  BackPredict(michel_pr). ToTH1(fdpot)->Write("FD_michel");
157  BackPredict(michelben_pr).ToTH1(fdpot)->Write("FD_michelben");
158  BackPredict(prop_pr). ToTH1(fdpot)->Write("FD_prop");
159 
160  SaveDCMPPlots(cheat, "MCEst", ndpot);
161  SaveDCMPPlots(&michel, "Michel", ndpot);
162  SaveDCMPPlots(&michelben, "MichelBen", ndpot);
163  SaveDCMPPlots(&prop, "Prop", ndpot);
164 
165 }
double GetRes(TH1 *h1, TH1 *h2, int j)
Near Detector underground.
Definition: SREnums.h:10
Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
virtual Spectrum AntiNueComponent() const =0
(&#39; appearance&#39;)
Definition: IPrediction.h:18
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:176
virtual Spectrum NumuComponent() const =0
HistAxis GetDCMPAxis(std::string PID)
virtual Spectrum NumuComponent() const override
osc::OscCalcDumb calc
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:39
const char * label
void SaveDCMPPlots(const IDecomp *dcmp, std::string label, double pot)
Charged-current interactions.
Definition: IPrediction.h:39
Uses MC for NC and CC components, assigns remainder of data to CC.
Definition: NumuDecomp.h:10
Cut GetPIDCut(std::string PID)
TFile * outFile
Definition: PlotXSec.C:135
Spectrum Predict(osc::IOscCalc *calc) const override
const Cut kNumuND
Definition: NumuCuts.h:55
const Var kCCE
Definition: NumuVars.h:21
#define pot
Simple oscillation probability calculator that has no solar term or mass hierarchy or delta so it&#39;s s...
Definition: OscCalcDumb.h:16
const double j
Definition: BetheBloch.cxx:29
virtual Spectrum AntiNumuComponent() const =0
double POT() const
Definition: Spectrum.h:226
bool CheckPIDString(std::string PID)
const SystShifts kNoShift
Definition: SystShifts.cxx:22
TH1F * h2
Definition: plot.C:45
TH1F * h1
void NueExtrap(string beam="fhc", string cvntype="oldpresel")
Definition: NueExtrap.C:28
code to link reconstructed objects back to the MC truth information
Definition: FillTruth.h:17
Splits Data proportionally according to MC.
const Var kTuftsWeightCC
Definition: XsecTunes.h:30
virtual Spectrum NCTotalComponent() const
Definition: IDecomp.h:18
std::vector< Loaders * > loaders
Definition: syst_header.h:386
Standard interface to all decomposition techniques.
Definition: IDecomp.h:13
void test_micheldecomp(std::string PID)
Spectrum BackPredict(PredictionExtrap pred)
assert(nhit_max >=nhit_nbins)
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
Take the output of an extrapolation and oscillate it as required.
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
Just return the ND truth spectra as the decomposition.
Definition: CheatDecomp.h:10
virtual Spectrum NueComponent() const =0
enum BeamMode string