11 #include "Utilities/func/MathUtil.h" 24 #include "TDirectory.h" 25 #include "TObjString.h" 48 :
fNC (NDMCloader, EAxis,
50 fAntiNC (NDMCloader, EAxis,
52 fNumu (NDMCloader, EAxis,
54 fAntiNumu (NDMCloader, EAxis,
56 fNue (NDMCloader, EAxis,
58 fAntiNue (NDMCloader, EAxis,
60 fData (NDDataloader, EAxis, cut, shiftData, wei),
61 fMichelNC (NDMCloader, EAxis,
MEAxis,
62 cut &&
kIsNC, shiftMC, wei),
63 fMichelNumu (NDMCloader, EAxis,
MEAxis,
65 fMichelNue (NDMCloader, EAxis,
MEAxis,
67 fMichelData (NDDataloader, EAxis,
MEAxis,
69 fNueEstimate(nueEstimate),
70 fOwnNueEstimate(false),
86 cut, nueEstimate, shiftMC, shiftData, wei)
98 axis, cut, nueEstimate, shiftMC, shiftData, wei)
108 std::vector<double> numuScales;
109 std::vector<double> ncScales;
110 std::vector<double> nueScales;
113 for (
int EIdx = 1; EIdx <= EIdxLimitHist->GetNbinsX(); EIdx++){
138 SplitByNue(hSums, numuScales, ncScales, nueScales);
142 delete EIdxLimitHist;
151 std::vector<double> &numuScales,
152 std::vector<double> &ncScales,
153 std::vector<double> &nueScales)
const 160 numuScales.push_back(numuScale);
161 ncScales.push_back(ncScale);
166 std::vector<double> &numuScales,
167 std::vector<double> &ncScales,
168 std::vector<double> &nueScales)
const 172 numuScales.push_back(1);
173 ncScales.push_back(1);
174 nueScales.push_back(1);
179 numuScales.push_back(DataMC);
180 ncScales.push_back(DataMC);
181 nueScales.push_back(DataMC);
185 std::vector<double> &numuScales,
186 std::vector<double> &ncScales,
187 std::vector<double> &nueScales)
const 199 TH1* hnuerat = nuerat.
ToTH1();
200 double scale = hnuerat->GetBinContent(EIdx);
210 TH1* hnumurat = numurat.
ToTH1();
211 double scale = hnumurat->GetBinContent(EIdx);
220 TH1* hncrat = ncrat.
ToTH1();
221 double scale = hncrat->GetBinContent(EIdx);
230 TH1* hnumurat = numurat.
ToTH1();
231 double numucont = hnumurat->GetBinContent(EIdx);
239 TH1* hncrat = ncrat.
ToTH1();
240 double nccont = hncrat->GetBinContent(EIdx);
255 double scale = (MaxScale-MinScale) / 2;
258 for (
double N = 0; N < 10; N++){
260 if ((MaxScale - MinScale) < 2*
dx)
261 dx = (MaxScale - MinScale)/4;
264 MaxScale = (likeplus > likeminus) ? MaxScale : scale;
265 MinScale = (likeplus > likeminus) ? scale : MinScale;
266 scale = (MaxScale+MinScale)/2;
273 assert(hSums.
fSumNumu > 0 &&
"Don't have enough stats. Bailing.");
285 assert(hSums.
fSumNC > 0 &&
"Don't have enough stats. Bailing.");
300 for (
int mIdx = 1; mIdx <= mIdxLimitHist->GetNbinsY(); mIdx++){
306 if (mean == 0 && Data == 0)
314 LL += Data - Data * TMath::Log(Data / mean) -
mean;
316 delete mIdxLimitHist;
324 double sum = hist->ProjectionX()->GetBinContent(EIdx);
331 int EIdx,
int mIdx)
const{
333 double content = hist->GetBinContent(EIdx, mIdx);
346 for (
int EIdx = 0; EIdx < scaledMC.size()-2; EIdx++){
349 double bincont = scales[EIdx]*scaledMC[binidx];
350 scaledMC[binidx] = bincont;
352 return Spectrum(std::move(scaledMC),
363 TH1* hdcmp = dcmp.
ToTH1(pot);
364 hdcmp->Write((comp+
"_micheldecomp").c_str());
366 hmc->Write((comp+
"_mc").c_str());
368 TH1* hratio = (TH1*)hdcmp->Clone(
"hratio");
370 hratio->Write((comp+
"_ratio").c_str());
372 delete hdcmp;
delete hmc;
delete hratio;
386 tempnumu->Write(
"NumuTemplate");
387 tempnc ->Write(
"NCTemplate");
388 tempnue ->Write(
"NueTemplate");
389 tempdata->Write(
"DataTemplate");
391 delete tempnumu;
delete tempnc;
delete tempnue;
delete tempdata;
397 TDirectory*
tmp = gDirectory;
508 TDirectory*
tmp = gDirectory;
510 dir = dir->mkdir(name.c_str());
513 TObjString(
"MichelDecomp").Write(
"type");
535 dir = dir->GetDirectory(name.c_str());
540 Spectrum* antinc = ana::LoadFrom<Spectrum>
544 Spectrum* antinumu = ana::LoadFrom<Spectrum>
548 Spectrum* antinue = ana::LoadFrom<Spectrum>
565 *numu, *antinumu, *nue, *antinue, *data,
566 *tempnc, *tempnumu, *tempnue, *tempdata, nueEst);
572 return std::unique_ptr<MichelDecomp>(
ret);
std::vector< double > ncScales
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
virtual Spectrum Data_Component() const override
_HistAxis< Var > HistAxis
double GetTemplateContent(ReweightableSpectrum spect, int EIdx, int mIdx) const
Gets bin content in template histogram, for a given E / NME bin.
Represent the binning of a Spectrum's x-axis.
Cuts and Vars for the 2020 FD DiF Study.
const Binning MEBinning
This is the Binning scheme used for the # Michel Templates.
Spectrum MC_AntiNumuComponent() const override
virtual Spectrum AntiNueComponent() const =0
TH1D * ToTH1(Color_t col=kBlack, Style_t style=kSolid) const
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.
virtual Spectrum NumuComponent() const =0
Simple record of shifts applied to systematic parameters.
double GetNumuScale(int EIdx) const
Numu scale from input decomp.
double GetNumuMCContent(int EIdx) const
Getting numu relative MC content.
double GetNCMCContent(int EIdx) const
Collection of SpectrumLoaders for many configurations.
virtual Spectrum AntiNumuComponent() const override
virtual Spectrum NCComponent() const override
Spectrum with the value of a second variable, allowing for reweighting
const Eigen::ArrayXd & GetEigen() const
NB these don't have POT scaling. For expert high performance ops only!
const Cut kIsAntiNu([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return sr->mc.nu[0].pdg< 0;})
Is this truly an antineutrino?
const Cut kIsBeamNue(CCFlavSel(12, 12))
Select CC .
virtual Spectrum MCToDCMPComp(std::vector< double >, Spectrum) const
Ending code to run through decomp algebra.
const IDecomp * fNueEstimate
double GetNCScale(int EIdx) const
NC scale from input decomp.
virtual Spectrum NumuComponent() const override
double MDCMPLogLikelihood(double h, int EIdx, MDCMPHelper hSums) const
Calculates the LL for a given scale, based on the data.
Representation of a spectrum in any variable, with associated POT.
virtual Spectrum AntiNueComponent() const override
static std::unique_ptr< MichelDecomp > LoadFrom(TDirectory *dir, const std::string &name)
void SavePlots(TDirectory *dir) const
const XML_Char const XML_Char * data
void Decompose() const
Calls ComputeScaleFactor for every process / E Bin, saves scale factors.
const std::vector< std::string > & GetLabels() const
void SaveTempPlots(TDirectory *dir) const
Spectrum MC_NCTotalComponent() const override
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
double CalculateMCMean(double h, double Numu, double NC, double Nue, MDCMPHelper hSums) const
Calculates the expected MC events, for a given scale factor.
ReweightableSpectrum fMichelNumu
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?
void SplitByMichels(int EIdx, MDCMPHelper hSums, std::vector< double > &numuScales, std::vector< double > &ncScales, std::vector< double > &nueScales) const
virtual Spectrum NueComponent() const override
ReweightableSpectrum fMichelData
Help reduce verbosity when passing hist values between MichelDecomp funcs.
void SaveTo(TDirectory *dir, const std::string &name) const
const std::vector< Binning > & GetBinnings() const
virtual Spectrum AntiNumuComponent() const =0
virtual Spectrum NCAntiComponent() const override
std::vector< float > Spectrum
std::vector< double > nueScales
std::vector< double > numuScales
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's assoc with slice.
double ComputeMaxScale(MDCMPHelper hSums) const
Represent the ratio between two spectra.
ReweightableSpectrum fMichelNue
Spectrum MC_NumuComponent() const override
Base class for the various types of spectrum loader.
Collect information describing the x-axis of an analysis histogram.
double ComputeScaleFactor(int EIdx, MDCMPHelper hSums) const
Starting code to run through decomp algebra /////////////////////////////////////////////////////////...
REGISTER_LOADFROM("BENDecomp", IDecomp, BENDecomp)
void SaveTo(TDirectory *dir, const std::string &name) const
virtual Spectrum NCTotalComponent() const
double GetSum(ReweightableSpectrum spect, int EIdx) const
Sums given spectrum over NME bins, for a given E Bin.
Standard interface to all decomposition techniques.
MichelDecomp(SpectrumLoaderBase &NDMCloader, SpectrumLoaderBase &NDdataloader, const HistAxis &axis, const Cut &cut, const IDecomp *nNueEstimate, const SystShifts &shiftMC=kNoShift, const SystShifts &shiftData=kNoShift, const Var &wei=kUnweighted)
assert(nhit_max >=nhit_nbins)
void SplitByProp(MDCMPHelper hSums, std::vector< double > &numuScales, std::vector< double > &ncScales, std::vector< double > &nueScales) const
Spectrum MC_NueComponent() const override
virtual void SaveTo(TDirectory *dir, const std::string &name) const =0
This module creates Common Analysis Files.
void SaveCompPlots(TDirectory *dir, std::string comp, Spectrum dcmp, Spectrum mc) const
const HistAxis MEAxis("N_{Michels}", MEBinning, kNMichels)
Spectrum MC_NCComponent() const override
std::unique_ptr< IDecomp > LoadFrom< IDecomp >(TDirectory *dir, const std::string &label)
MDCMPFitResults fFitResult
Spectrum MC_AntiNueComponent() const override
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
double GetNueScale(int EIdx) const
Nue scale from input decomp.
TH2 * ToTH2(const Spectrum &s, double exposure, ana::EExposureType expotype, const Binning &binsx, const Binning &binsy, ana::EBinType bintype)
For use with Var2D.
virtual Spectrum NueComponent() const =0
ReweightableSpectrum fMichelNC
virtual Spectrum NCTotalComponent() const override
TH2D * ToTH2(double pot) const
void SplitByNue(MDCMPHelper hSums, std::vector< double > &numuScales, std::vector< double > &ncScales, std::vector< double > &nueScales) const
Spectrum MC_NCAntiComponent() const override