MichelDecompTest.C
Go to the documentation of this file.
2 
3 #include "CAFAna/Cuts/Cuts.h"
5 #include "CAFAna/Core/Spectrum.h"
11 #include "CAFAna/Analysis/Plots.h"
12 #include "CAFAna/Systs/Systs.h"
13 #include "CAFAna/Core/SystShifts.h"
15 #include "CAFAna/Vars/Vars.h"
16 #include "CAFAna/Fit/Fit.h"
17 #include "CAFAna/Vars/FitVars.h"
18 #include "CAFAna/Core/Loaders.h"
22 
23 #include "Utilities/func/MathUtil.h"
24 
25 #include "TCanvas.h"
26 #include "TH2.h"
27 #include "TFile.h"
28 #include "TLegend.h"
29 
30 #include <iostream>
31 #include <cmath>
32 
33 using namespace ana;
34 
36 {
37  calc->SetL(810);
38  calc->SetRho(2.75);
39  calc->SetDmsq21(7.59e-5);
40  calc->SetDmsq32(2.40e-3);
41  calc->SetTh12(asin(sqrt(.87))/2);
42  calc->SetTh13(asin(sqrt(.095))/2);
43  calc->SetTh23(asin(sqrt(.95))/2);
44  calc->SetdCP(0);
45 }
46 
48 {
49 
50  const Cut cutNDContNue( kNDContainNue );
51  const Cut cutNDPIDNue( kLEM > .64 );
52  const Cut cutNDNue( cutNDContNue && cutNDPIDNue );
53 
54  const Binning MBinning = Binning::Simple(8,0,8);
55 
56  const std::string farNonS = "$NOVA_PROD/mc/FA14-10-28/genie/fd/caf/000148/*/*fhc_nonswap*.caf.root";
57  const std::string farSwap = "$NOVA_PROD/mc/FA14-10-28/genie/fd/caf/000148/*/*fhc_swap*.caf.root";
58  const std::string nearNonS = "$NOVA_PROD/mc/FA14-11-25/genie/nd/caf/*103/10378/*s101*.caf.root";
59 
60  SpectrumLoader loaderNonS(farNonS);
61  SpectrumLoader loaderSwap(farSwap);
62  SpectrumLoader loaderNDMC(nearNonS);
63 
64  const Binning bins = Binning::Simple(6, 1.25, 2.75);
65  const HistAxis EAxis("Energy (GeV)", bins, kCaloE);
66 
67  SystShifts shift;
68  shift.SetShift(new GenieRwgtTableSyst(rwgt::fReweightRvnCC2pi), +2);
69  shift.SetShift(new GenieRwgtTableSyst(rwgt::fReweightRvnNC2pi), +2);
70  shift.SetShift(new GenieRwgtTableSyst(rwgt::fReweightRvbarnCC2pi), +2);
71  shift.SetShift(new GenieRwgtTableSyst(rwgt::fReweightRvbarnNC2pi), +2);
72 
73  Spectrum CCMC (loaderNDMC, EAxis, cutNDNue && kIsNumuCC);
74  Spectrum NCMC (loaderNDMC, EAxis, cutNDNue && kIsNC);
75  Spectrum NEMC (loaderNDMC, EAxis, cutNDNue && kIsBeamNue);
76 
77  Spectrum CCData(loaderNDMC, EAxis, cutNDNue && kIsNumuCC, shift);
78  Spectrum NCData(loaderNDMC, EAxis, cutNDNue && kIsNC, shift);
79  Spectrum NEData(loaderNDMC, EAxis, cutNDNue && kIsBeamNue, shift);
80 
81  // N Michels kNoShifted / nominal spectra
82  Spectrum CCMCME (EAxis.label, loaderNDMC,
83  EAxis.bins, EAxis.var,
84  MBinning, kNMichels,
85  cutNDNue && kIsNumuCC);
86  Spectrum NCMCME (EAxis.label, loaderNDMC,
87  EAxis.bins, EAxis.var,
88  MBinning, kNMichels,
89  cutNDNue && kIsNC);
90  Spectrum NEMCME (EAxis.label, loaderNDMC,
91  EAxis.bins, EAxis.var,
92  MBinning, kNMichels,
93  cutNDNue && kIsBeamNue);
94  Spectrum CCDataME(EAxis.label, loaderNDMC,
95  EAxis.bins, EAxis.var,
96  MBinning, kNMichels,
97  cutNDNue && kIsNumuCC, shift);
98  Spectrum NCDataME(EAxis.label, loaderNDMC,
99  EAxis.bins, EAxis.var,
100  MBinning, kNMichels,
101  cutNDNue && kIsNC, shift);
102  Spectrum NEDataME(EAxis.label, loaderNDMC,
103  EAxis.bins, EAxis.var,
104  MBinning, kNMichels,
105  cutNDNue && kIsBeamNue, shift);
106 
107  MichelDecomp mDCMP(loaderNDMC, loaderNDMC,
108  EAxis, cutNDNue, kNoShift, shift);
109  NumuDecomp numuDCMP(loaderNDMC, loaderNDMC,
110  EAxis, kNDContainNue && kNumuCC, kNoShift, shift);
111  NueExtrap extrapDCMP(loaderNDMC,
112  loaderSwap,
113  loaderNonS,
114  kNullLoader,
115  mDCMP,
116  numuDCMP,
117  EAxis,
118  EAxis,
119  cutNDPIDNue,
120  cutNDNue,
122  kNoShift);
123  PredictionNoExtrap predNominal(loaderNonS, loaderSwap,
124  EAxis.label, EAxis.bins, EAxis.var, cutNDPIDNue, kNoShift);
125 
126  PredictionExtrap predDCMP(&extrapDCMP);
127 
128  PredictionNoExtrap predShifts(loaderNonS, loaderSwap,
129  EAxis.label, EAxis.bins, EAxis.var, cutNDPIDNue, shift);
130 
131 
132  // Decomp Loaders
133  loaderNDMC.Go();
134  loaderSwap.Go();
135  loaderNonS.Go();
136 
137  Spectrum CCDCMP = mDCMP.NumuComponent() + mDCMP.AntiNumuComponent();
138  Spectrum NCDCMP = mDCMP.NCTotalComponent();
139  Spectrum NEDCMP = mDCMP.NueComponent() + mDCMP.AntiNueComponent();
140 
142  resetCalc(&calc);
143 
144  double ourPOT = CCDataME.POT();
145 
146  Spectrum obsNominal = predNominal.Predict(&calc).FakeData(ourPOT);
147  SingleSampleExperiment exptNominal(&predNominal, obsNominal);
148 
149  Spectrum obsShifts = predShifts.Predict(&calc).FakeData(ourPOT);
150  SingleSampleExperiment exptShifts(&predNominal, obsShifts);
151 
152  Spectrum obsDCMP = predShifts.Predict(&calc).FakeData(ourPOT);
153  SingleSampleExperiment exptDCMP(&predDCMP, obsDCMP);
154 
155  Spectrum obsBlah = predDCMP.Predict(&calc).FakeData(ourPOT);
156 
157  Spectrum bkgdNominal = obsNominal - predNominal.PredictComponent(&calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kBoth);
158  Spectrum bkgdShifted = obsShifts - predShifts.PredictComponent(&calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kBoth);
159  Spectrum bkgdDCMP = obsBlah - predDCMP.PredictComponent(&calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kBoth);
160 
161  TFile* outFile = new TFile("Rvn2pi.root","RECREATE");
162  outFile->cd();
163 
164  CCMCME.ToTH2(ourPOT)->Write("ME_Dist_CCMC");
165  NCMCME.ToTH2(ourPOT)->Write("ME_Dist_NCMC");
166  NEMCME.ToTH2(ourPOT)->Write("ME_Dist_NuEMC");
167  CCDataME.ToTH2(ourPOT)->Write("ME_Dist_CCData");
168  NCDataME.ToTH2(ourPOT)->Write("ME_Dist_NCData");
169  NEDataME.ToTH2(ourPOT)->Write("ME_Dist_NuEData");
170 
171  obsBlah.ToTH1(ourPOT)->Write("obsDCMP");
172  obsNominal.ToTH1(ourPOT)->Write("obsNominal");
173  obsShifts.ToTH1(ourPOT)->Write("obsShifted");
174 
175  bkgdDCMP.ToTH1(ourPOT)->Write("bkgdDCMP");
176  bkgdNominal.ToTH1(ourPOT)->Write("bkgdNominal");
177  bkgdShifted.ToTH1(ourPOT)->Write("bkgdShifted");
178 
179  CCMC.ToTH1(ourPOT)->Write("ND_CCMC");
180  CCData.ToTH1(ourPOT)->Write("ND_CCData");
181  CCDCMP.ToTH1(ourPOT)->Write("ND_CCDCMP");
182 
183  NCMC.ToTH1(ourPOT)->Write("ND_NCMC");
184  NCData.ToTH1(ourPOT)->Write("ND_NCData");
185  NCDCMP.ToTH1(ourPOT)->Write("ND_NCDCMP");
186 
187  NEMC.ToTH1(ourPOT)->Write("ND_NEMC");
188  NEData.ToTH1(ourPOT)->Write("ND_NEData");
189  NEDCMP.ToTH1(ourPOT)->Write("ND_NEDCMP");
190 
191 }
Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
const Var kLEM
PID
Definition: Vars.cxx:26
virtual void SetL(double L)=0
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
Definition: Spectrum.cxx:193
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void MichelDecompTest()
(&#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 void SetDmsq21(const T &dmsq21)=0
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
T sqrt(T number)
Definition: d0nt_math.hpp:156
virtual Spectrum AntiNumuComponent() const override
const Cut kNumuCC
Definition: NumuCuts.h:70
virtual void SetTh13(const T &th13)=0
const Cut kIsBeamNue(CCFlavSel(12, 12))
Select CC .
virtual Spectrum NumuComponent() const override
osc::OscCalcDumb calc
virtual void SetDmsq32(const T &dmsq32)=0
void resetCalc(osc::IOscCalcAdjustable *calc)
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:39
virtual Spectrum AntiNueComponent() const override
Charged-current interactions.
Definition: IPrediction.h:39
Uses MC for NC and CC components, assigns remainder of data to CC.
Definition: NumuDecomp.h:10
TFile * outFile
Definition: PlotXSec.C:135
const Var kCaloE
Summed calorimetric energy of all hits in slice, uncorrected.
Definition: Vars.cxx:52
Spectrum Predict(osc::IOscCalc *calc) const override
const Cut kIsNC([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return!sr->mc.nu[0].iscc;})
Is this a Neutral Current event?
Definition: TruthCuts.h:8
Spectrum FakeData(double pot) const
Synonymous with AsimovData(). Retained for compatibility.
Definition: Spectrum.cxx:377
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:32
virtual Spectrum NueComponent() const override
virtual void Go() override
Load all the registered spectra.
double POT() const
Definition: Spectrum.h:226
const Var kNMichels([](const caf::SRProxy *sr){int n_me=0;for(int i=0;i< (int) sr->me.nslc;i++) if(sr->me.slc[i].mid > 1. &&sr->me.slc[i].deltat > 1200.) n_me++;for(int i=0;i< (int) sr->me.nkalman;i++) if(sr->me.trkkalman[i].mid > 0. &&sr->me.trkkalman[i].deltat > 800.) n_me++;if(n_me > 2) n_me=2;return n_me;})
Number of SlcME&#39;s assoc with slice.
Definition: Vars.h:85
const SystShifts kNoShift
Definition: SystShifts.cxx:22
const Binning bins
Definition: NumuCC_CPiBin.h:8
const Cut kNDContainNue([](const caf::SRProxy *sr){if(sr->vtx.elastic.fuzzyk.nshwlid<=0) return false;const caf::SRShowerLIDProxy &shwlid=sr->vtx.elastic.fuzzyk.png[0].shwlid;const caf::SRVector3DProxy &start=shwlid.start;const caf::SRVector3DProxy &stop=shwlid.stop;return(std::min(start.X(), stop.X()) >-150 && std::max(start.X(), stop.X())< +150 && std::min(start.Y(), stop.Y()) >-150 && std::max(start.Y(), stop.Y())< +150 && std::min(start.Z(), stop.Z()) > 50 && std::max(start.Z(), stop.Z())< 1230);})
Definition: Cuts.h:22
virtual void SetRho(double rho)=0
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
virtual void SetTh23(const T &th23)=0
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Prediction that just uses FD MC, with no extrapolation.
Take the output of an extrapolation and oscillate it as required.
Float_t e
Definition: plot.C:35
loaderSwap
Definition: demo4.py:9
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:81
static NullLoader kNullLoader
Dummy loader that doesn&#39;t load any files.
virtual void SetTh12(const T &th12)=0
virtual Spectrum NCTotalComponent() const override
virtual void SetdCP(const T &dCP)=0
T asin(T number)
Definition: d0nt_math.hpp:60
Compare a single data spectrum to the MC + cosmics expectation.
enum BeamMode string