MichelDecomp.cxx
Go to the documentation of this file.
3 #include "CAFAna/Core/Ratio.h"
7 #include "CAFAna/Core/Loaders.h"
9 #include "CAFAna/Vars/Vars.h"
10 
11 #include "Utilities/func/MathUtil.h"
12 
13 #include <cassert>
14 
15 #include "TH1.h"
16 #include "TH2.h"
17 #include "TH3.h"
18 #include "TCanvas.h"
19 #include "TProfile.h"
20 #include "TAxis.h"
21 #include "TMatrixD.h"
22 #include "TMath.h"
23 #include "TVectorD.h"
24 #include "TDirectory.h"
25 #include "TObjString.h"
26 #include "TRandom.h"
27 
28 namespace ana
29 {
30  REGISTER_LOADFROM("MichelDecomp", IDecomp, MichelDecomp);
31 
32  //----------------------------------------------------------------------
33 
34 
36  const HistAxis MEAxis("N_{Michels}",MEBinning,kNMichels);
37 
38  //----------------------------------------------------------------------
39  //----------------------------------------------------------------------
41  SpectrumLoaderBase& NDDataloader,
42  const HistAxis& EAxis,
43  const Cut& cut,
44  const IDecomp* nueEstimate,
45  const SystShifts& shiftMC,
46  const SystShifts& shiftData,
47  const Var& wei)
48  : fNC (NDMCloader, EAxis,
49  cut && kIsNC && !kIsAntiNu, shiftMC, wei),
50  fAntiNC (NDMCloader, EAxis,
51  cut && kIsNC && kIsAntiNu, shiftMC, wei),
52  fNumu (NDMCloader, EAxis,
53  cut && kIsNumuCC && !kIsAntiNu, shiftMC, wei),
54  fAntiNumu (NDMCloader, EAxis,
55  cut && kIsNumuCC && kIsAntiNu, shiftMC, wei),
56  fNue (NDMCloader, EAxis,
57  cut && kIsBeamNue && !kIsAntiNu, shiftMC, wei),
58  fAntiNue (NDMCloader, EAxis,
59  cut && kIsBeamNue && kIsAntiNu, shiftMC, wei),
60  fData (NDDataloader, EAxis, cut, shiftData, wei),
61  fMichelNC (NDMCloader, EAxis, MEAxis,
62  cut && kIsNC, shiftMC, wei),
63  fMichelNumu (NDMCloader, EAxis, MEAxis,
64  cut && kIsNumuCC, shiftMC, wei),
65  fMichelNue (NDMCloader, EAxis, MEAxis,
66  cut && kIsBeamNue, shiftMC, wei),
67  fMichelData (NDDataloader, EAxis, MEAxis,
68  cut, shiftData, wei),
69  fNueEstimate(nueEstimate),
70  fOwnNueEstimate(false),
71  isDecomposed(false)
72  {}
73 
75  SpectrumLoaderBase& NDDataloader,
77  const Binning& bins,
78  const Var& var,
79  const Cut& cut,
80  const IDecomp* nueEstimate,
81  const SystShifts& shiftMC,
82  const SystShifts& shiftData,
83  const Var& wei)
84  : MichelDecomp(NDMCloader, NDDataloader,
85  HistAxis(label, bins, var),
86  cut, nueEstimate, shiftMC, shiftData, wei)
87  {}
88 
90  const HistAxis& axis,
91  const Cut& cut,
92  const IDecomp* nueEstimate,
93  const SystShifts& shiftMC,
94  const SystShifts& shiftData,
95  const Var& wei)
96  : MichelDecomp(loaders.GetLoader(caf::kNEARDET, Loaders::kMC),
97  loaders.GetLoader(caf::kNEARDET, Loaders::kData),
98  axis, cut, nueEstimate, shiftMC, shiftData, wei)
99  {}
100 
102  {
103  if(fOwnNueEstimate) delete fNueEstimate;
104  }
105 
106  /// Calls ComputeScaleFactor for every process / E Bin, saves scale factors
108  std::vector<double> numuScales;
109  std::vector<double> ncScales;
110  std::vector<double> nueScales;
111 
112  TH1* EIdxLimitHist = fNumu.ToTH1(1e20);
113  for (int EIdx = 1; EIdx <= EIdxLimitHist->GetNbinsX(); EIdx++){
114  // Get appropriate nuescale from input decomp
115 
116  // Stores sums of our histograms, to reduce computing in functions
117  MDCMPHelper hSums;
118  hSums.fSumData = GetSum(fMichelData, EIdx);
119  hSums.fSumNumu = GetSum(fMichelNumu, EIdx);
120  hSums.fSumNC = GetSum(fMichelNC, EIdx);
121  hSums.fSumNue = GetSum(fMichelNue, EIdx);
122  hSums.fNueScale = GetNueScale(EIdx);
123  hSums.fNumuScale = GetNumuScale(EIdx);
124  hSums.fNCScale = GetNCScale(EIdx);
125 
126  // Typical nue slices have about 1/8 SlcME's per slice
127  // So, if there are more than about 80 events in a bin, we have
128  // at least 10 SlcME's in that bin, which is enough to run MDCMP
129  // If we have fewer, proportionally split NC's and CC's
130 
131  // Instead of cheacking total number of events in a bin, look at
132  // number of data events seen, less the nue component predicted
133  /*if (EIdx == 1 || EIdx == 2 || EIdx == 8 || EIdx == 9 || EIdx == 10 ||
134  EIdx ==11 || EIdx ==16 || EIdx ==17 || EIdx ==18 || EIdx == 19 ||
135  EIdx ==20 || EIdx ==21 || EIdx > 23)*/
136  if ((GetNumuMCContent(EIdx) < .2) || (GetNCMCContent(EIdx) == .0))
137  //SplitByProp(hSums, numuScales, ncScales, nueScales);
138  SplitByNue(hSums, numuScales, ncScales, nueScales);
139  else
140  SplitByMichels(EIdx, hSums, numuScales, ncScales, nueScales);
141  }
142  delete EIdxLimitHist;
143 
144  fFitResult.numuScales = numuScales;
145  fFitResult.ncScales = ncScales;
146  fFitResult.nueScales = nueScales;
147  isDecomposed = true;
148  }
149 
151  std::vector<double> &numuScales,
152  std::vector<double> &ncScales,
153  std::vector<double> &nueScales) const
154  {
155  // We have enough data in this bin to do MichelDecomp
156  double numuScale = ComputeScaleFactor(EIdx, hSums);
157  double ncScale =
158  (hSums.fSumData - numuScale*hSums.fSumNumu -
159  hSums.fNueScale*hSums.fSumNue) / hSums.fSumNC;
160  numuScales.push_back(numuScale);
161  ncScales.push_back(ncScale);
162  nueScales.push_back(hSums.fNueScale);
163  }
164 
166  std::vector<double> &numuScales,
167  std::vector<double> &ncScales,
168  std::vector<double> &nueScales) const
169  {
170  // Return MC if there are no events in the MC
171  if (hSums.fSumNumu + hSums.fSumNC + hSums.fSumNue == 0){
172  numuScales.push_back(1);
173  ncScales.push_back(1);
174  nueScales.push_back(1);
175  return;
176  }
177  double DataMC = hSums.fSumData /
178  (hSums.fSumNue + hSums.fSumNumu + hSums.fSumNC);
179  numuScales.push_back(DataMC);
180  ncScales.push_back(DataMC);
181  nueScales.push_back(DataMC);
182  }
183 
185  std::vector<double> &numuScales,
186  std::vector<double> &ncScales,
187  std::vector<double> &nueScales) const
188  {
189  numuScales.push_back(hSums.fNumuScale);
190  ncScales.push_back(hSums.fNCScale);
191  nueScales.push_back(hSums.fNueScale);
192  }
193 
194  /// Nue scale from input decomp
195  double MichelDecomp::GetNueScale(int EIdx) const {
198  Ratio nuerat = nue / (fNue+fAntiNue);
199  TH1* hnuerat = nuerat.ToTH1();
200  double scale = hnuerat->GetBinContent(EIdx);
201  delete hnuerat;
202  return scale;
203  }
204 
205  /// Numu scale from input decomp
206  double MichelDecomp::GetNumuScale(int EIdx) const {
209  Ratio numurat = numu / (fNumu+fAntiNumu);
210  TH1* hnumurat = numurat.ToTH1();
211  double scale = hnumurat->GetBinContent(EIdx);
212  delete hnumurat;
213  return scale;
214  }
215 
216  /// NC scale from input decomp
217  double MichelDecomp::GetNCScale(int EIdx) const {
219  Ratio ncrat = nc / (fNC+fAntiNC);
220  TH1* hncrat = ncrat.ToTH1();
221  double scale = hncrat->GetBinContent(EIdx);
222  delete hncrat;
223  return scale;
224  }
225 
226  /// Getting numu relative MC content
227  double MichelDecomp::GetNumuMCContent(int EIdx) const {
229  Ratio numurat = numu / (numu+fNC+fAntiNC+fNue+fAntiNue);
230  TH1* hnumurat = numurat.ToTH1();
231  double numucont = hnumurat->GetBinContent(EIdx);
232  delete hnumurat;
233  return numucont;
234  }
235 
236  double MichelDecomp::GetNCMCContent(int EIdx) const {
238  Ratio ncrat = nc / (nc+fNumu+fAntiNumu+fNue+fAntiNue);
239  TH1* hncrat = ncrat.ToTH1();
240  double nccont = hncrat->GetBinContent(EIdx);
241  delete hncrat;
242  return nccont;
243  }
244 
245  ///////////////////////////////////////////////////////////////////////////
246  /// Starting code to run through decomp algebra
247  ///////////////////////////////////////////////////////////////////////////
248  /// Compute scale factors
249  double MichelDecomp::ComputeScaleFactor(int EIdx, MDCMPHelper hSums) const {
250  double dx = 0.01;
251 
252  // Gives the bounds for the scale factor after each iteration
253  double MinScale = 0;
254  double MaxScale = ComputeMaxScale(hSums);
255  double scale = (MaxScale-MinScale) / 2;
256 
257  // Continue for 10 iterations
258  for (double N = 0; N < 10; N++){
259  // Adjust dx to be within the current allowed scale factors
260  if ((MaxScale - MinScale) < 2*dx)
261  dx = (MaxScale - MinScale)/4;
262  double likeplus = MDCMPLogLikelihood(scale+dx, EIdx, hSums);
263  double likeminus = MDCMPLogLikelihood(scale-dx, EIdx, hSums);
264  MaxScale = (likeplus > likeminus) ? MaxScale : scale;
265  MinScale = (likeplus > likeminus) ? scale : MinScale;
266  scale = (MaxScale+MinScale)/2;
267  }
268  return scale;
269  }
270 
271 
273  assert(hSums.fSumNumu > 0 && "Don't have enough stats. Bailing.");
274  double MaxScale = (hSums.fSumData - hSums.fNueScale*hSums.fSumNue) /
275  hSums.fSumNumu;
276  return MaxScale;
277  }
278 
279  /// Calculates the expected MC events, for a given scale factor
281  double Numu,
282  double NC,
283  double Nue,
284  MDCMPHelper hSums) const{
285  assert(hSums.fSumNC > 0 && "Don't have enough stats. Bailing.");
286 
287  double mean =
288  h*(Numu - hSums.fSumNumu / hSums.fSumNC * NC) +
289  hSums.fSumData / hSums.fSumNC * NC -
290  hSums.fNueScale*hSums.fSumNue / hSums.fSumNC * NC + hSums.fNueScale*Nue;
291  return mean;
292  }
293 
294  /// Calculates the LL for a given scale, based on the data
295  double MichelDecomp::MDCMPLogLikelihood(double h, int EIdx,
296  MDCMPHelper hSums) const{
297  double LL = 0;
298 
299  TH1* mIdxLimitHist = fMichelNumu.ToTH2(1e20);
300  for (int mIdx = 1; mIdx <= mIdxLimitHist->GetNbinsY(); mIdx++){
301  double Data = GetTemplateContent(fMichelData, EIdx, mIdx);
302  double Numu = GetTemplateContent(fMichelNumu, EIdx, mIdx);
303  double NC = GetTemplateContent(fMichelNC, EIdx, mIdx);
304  double Nue = GetTemplateContent(fMichelNue, EIdx, mIdx);
305  double mean = CalculateMCMean(h, Numu, NC, Nue, hSums);
306  if (mean == 0 && Data == 0)
307  continue;
308  if (Data == 0){
309  LL += Data - mean;
310  continue;
311  }
312  if (mean == 0)
313  mean = 1e-10;
314  LL += Data - Data * TMath::Log(Data / mean) - mean;
315  }
316  delete mIdxLimitHist;
317 
318  return LL;
319  }
320 
321  /// Sums given spectrum over NME bins, for a given E Bin
322  double MichelDecomp::GetSum(ReweightableSpectrum spect, int EIdx) const{
323  TH2* hist = spect.ToTH2(fData.POT());
324  double sum = hist->ProjectionX()->GetBinContent(EIdx);
325  delete hist;
326  return sum;
327  }
328 
329  /// Gets bin content in template histogram, for a given E / NME bin.
331  int EIdx, int mIdx) const{
332  TH2* hist = spec.ToTH2(fData.POT());
333  double content = hist->GetBinContent(EIdx, mIdx);
334  delete hist;
335  return content;
336  }
337  ///////////////////////////////////////////////////////////////////////////
338  /// Ending code to run through decomp algebra
339  ///////////////////////////////////////////////////////////////////////////
340 
341  // Inputs a std::vector of scale factors, and a Comp spectrum,
342  // And outputs decomp result with components scaled by MDCMP fit
343  Spectrum MichelDecomp::MCToDCMPComp(std::vector<double> scales,
344  Spectrum MCComp) const{
345  Eigen::ArrayXd scaledMC = MCComp.GetEigen(fData.POT());
346  for (int EIdx = 0; EIdx < scaledMC.size()-2; EIdx++){
347  // Bin content in new histogram is old content * fit scale
348  int binidx = EIdx+1; // Root counts bins starting at 1
349  double bincont = scales[EIdx]*scaledMC[binidx];
350  scaledMC[binidx] = bincont;
351  }
352  return Spectrum(std::move(scaledMC),
354  fData.POT(), fData.Livetime());
355  }
356 
357  //----------------------------------------------------------------------
359  Spectrum dcmp, Spectrum mc) const
360  {
361  double pot = fData.POT();
362  dir->cd();
363  TH1* hdcmp = dcmp.ToTH1(pot);
364  hdcmp->Write((comp+"_micheldecomp").c_str());
365  TH1* hmc = mc.ToTH1(pot);
366  hmc->Write((comp+"_mc").c_str());
367 
368  TH1* hratio = (TH1*)hdcmp->Clone("hratio");
369  hratio->Divide(hmc);
370  hratio->Write((comp+"_ratio").c_str());
371 
372  delete hdcmp; delete hmc; delete hratio;
373  }
374 
375  //----------------------------------------------------------------------
376  void MichelDecomp::SaveTempPlots(TDirectory* dir) const
377  {
378  double pot = fData.POT();
379  dir->cd();
380 
381  TH2* tempnumu = fMichelNumu.ToTH2(pot);
382  TH2* tempnc = fMichelNC. ToTH2(pot);
383  TH2* tempnue = fMichelNue. ToTH2(pot);
384  TH2* tempdata = fMichelData.ToTH2(pot);
385 
386  tempnumu->Write("NumuTemplate");
387  tempnc ->Write("NCTemplate");
388  tempnue ->Write("NueTemplate");
389  tempdata->Write("DataTemplate");
390 
391  delete tempnumu; delete tempnc; delete tempnue; delete tempdata;
392  }
393 
394  //----------------------------------------------------------------------
395  void MichelDecomp::SavePlots(TDirectory* dir) const
396  {
397  TDirectory* tmp = gDirectory;
398  dir->cd();
399 
400  SaveCompPlots(dir->mkdir("numu"), "numu",
401  NumuComponent(), fNumu);
402  SaveCompPlots(dir->mkdir("antinumu"), "antinumu",
404  SaveCompPlots(dir->mkdir("nc"), "nc",
405  NCComponent(), fNC);
406  SaveCompPlots(dir->mkdir("antinc"), "antinc",
408  SaveCompPlots(dir->mkdir("nue"), "nue",
409  NueComponent(), fNue);
410  SaveCompPlots(dir->mkdir("antinue"), "antinue",
412 
413  SaveTempPlots(dir->mkdir("templates"));
414 
415  tmp->cd();
416  }
417 
418  //----------------------------------------------------------------------
420  if (!isDecomposed)
421  Decompose();
423  }
424 
425  //----------------------------------------------------------------------
427  if (!isDecomposed)
428  Decompose();
430  }
431 
432  //----------------------------------------------------------------------
434  if (!isDecomposed)
435  Decompose();
437  }
438 
439  //----------------------------------------------------------------------
441  if (!isDecomposed)
442  Decompose();
444  }
445 
446  //----------------------------------------------------------------------
448  if (!isDecomposed)
449  Decompose();
450  return NCComponent()+NCAntiComponent();
451  }
452 
454  if (!isDecomposed)
455  Decompose();
457  }
458 
460  if (!isDecomposed)
461  Decompose();
463  }
464 
465  // And now make the MC distributions publicly accessible
466  //----------------------------------------------------------------------
468  return fNue;
469  }
470 
471  //----------------------------------------------------------------------
473  return fAntiNue;
474  }
475 
476  //----------------------------------------------------------------------
478  return fNumu;
479  }
480 
481  //----------------------------------------------------------------------
483  return fAntiNumu;
484  }
485 
486  //----------------------------------------------------------------------
488  return fNC+fAntiNC;
489  }
490 
491  //----------------------------------------------------------------------
493  return fNC;
494  }
495 
496  //----------------------------------------------------------------------
498  return fAntiNC;
499  }
500 
501  //----------------------------------------------------------------------
503  return fData;
504  }
505 
506  //----------------------------------------------------------------------
507  void MichelDecomp::SaveTo(TDirectory* dir, const std::string& name) const{
508  TDirectory* tmp = gDirectory;
509 
510  dir = dir->mkdir(name.c_str()); // switch to subdir
511  dir->cd();
512 
513  TObjString("MichelDecomp").Write("type");
514 
515  fNC.SaveTo( dir, "NC");
516  fAntiNC.SaveTo( dir, "AntiNC");
517  fNumu.SaveTo( dir, "Numu");
518  fAntiNumu.SaveTo( dir, "AntiNumu");
519  fNue.SaveTo( dir, "Nue");
520  fAntiNue.SaveTo( dir, "AntiNue");
521  fData.SaveTo( dir, "Data");
522  fMichelNC.SaveTo( dir, "TempNC");
523  fMichelNumu.SaveTo( dir, "TempNumu");
524  fMichelNue.SaveTo( dir, "TempNue");
525  fMichelData.SaveTo( dir, "TempData");
526  fNueEstimate->SaveTo(dir, "NueEstimate");
527 
528  dir->Write();
529  delete dir;
530 
531  tmp->cd();
532  }
533 
534  std::unique_ptr<MichelDecomp> MichelDecomp::LoadFrom(TDirectory* dir, const std::string& name){
535  dir = dir->GetDirectory(name.c_str()); // switch to subdir
536  assert(dir);
537 
538  Spectrum* nc = ana::LoadFrom<Spectrum>
539  (dir, "NC").release();
540  Spectrum* antinc = ana::LoadFrom<Spectrum>
541  (dir, "AntiNC").release();
542  Spectrum* numu = ana::LoadFrom<Spectrum>
543  (dir, "Numu").release();
544  Spectrum* antinumu = ana::LoadFrom<Spectrum>
545  (dir, "AntiNumu").release();
546  Spectrum* nue = ana::LoadFrom<Spectrum>
547  (dir, "Nue").release();
548  Spectrum* antinue = ana::LoadFrom<Spectrum>
549  (dir, "AntiNue").release();
550  Spectrum* data = ana::LoadFrom<Spectrum>
551  (dir, "Data").release();
552  ReweightableSpectrum* tempnc = ana::LoadFrom<ReweightableSpectrum>
553  (dir, "TempNC").release();
554  ReweightableSpectrum* tempnumu = ana::LoadFrom<ReweightableSpectrum>
555  (dir, "TempNumu").release();
556  ReweightableSpectrum* tempnue = ana::LoadFrom<ReweightableSpectrum>
557  (dir, "TempNue").release();
558  ReweightableSpectrum* tempdata = ana::LoadFrom<ReweightableSpectrum>
559  (dir, "TempData").release();
560  const IDecomp* nueEst = ana::LoadFrom<IDecomp>
561  (dir, "NueEstimate").release();
562 
563  MichelDecomp* ret =
564  new MichelDecomp(*nc, *antinc,
565  *numu, *antinumu, *nue, *antinue, *data,
566  *tempnc, *tempnumu, *tempnue, *tempdata, nueEst);
567 
568  ret->fOwnNueEstimate = true;
569 
570  delete dir;
571 
572  return std::unique_ptr<MichelDecomp>(ret);
573  }
574 }
575 
576 // LocalWords: EIdx
std::vector< double > ncScales
Definition: MichelDecomp.h:26
const XML_Char * name
Definition: expat.h:151
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
virtual Spectrum Data_Component() const override
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
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&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
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
Definition: Ratio.cxx:67
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:148
virtual Spectrum NumuComponent() const =0
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
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.
Definition: Loaders.h:23
virtual Spectrum AntiNumuComponent() const override
virtual Spectrum NCComponent() const override
const Color_t kMC
Spectrum with the value of a second variable, allowing for reweighting
const Eigen::ArrayXd & GetEigen() const
NB these don&#39;t have POT scaling. For expert high performance ops only!
Definition: Spectrum.h:188
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?
Definition: TruthCuts.h:53
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
Definition: MichelDecomp.h:195
Float_t tmp
Definition: plot.C:36
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.
Definition: Spectrum.h:40
virtual Spectrum AntiNueComponent() const override
static std::unique_ptr< MichelDecomp > LoadFrom(TDirectory *dir, const std::string &name)
void SavePlots(TDirectory *dir) const
const char * label
const XML_Char const XML_Char * data
Definition: expat.h:268
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
Double_t scale
Definition: plot.C:25
Spectrum MC_NCTotalComponent() const override
double dx[NP][NC]
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
Definition: MichelDecomp.h:191
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
void SplitByMichels(int EIdx, MDCMPHelper hSums, std::vector< double > &numuScales, std::vector< double > &ncScales, std::vector< double > &nueScales) const
virtual Spectrum NueComponent() const override
#define pot
ReweightableSpectrum fMichelData
Definition: MichelDecomp.h:193
Help reduce verbosity when passing hist values between MichelDecomp funcs.
Definition: MichelDecomp.h:31
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:517
const std::vector< Binning > & GetBinnings() const
virtual Spectrum AntiNumuComponent() const =0
virtual Spectrum NCAntiComponent() const override
std::vector< float > Spectrum
Definition: Constants.h:570
std::vector< double > nueScales
Definition: MichelDecomp.h:27
std::vector< double > numuScales
Definition: MichelDecomp.h:25
double POT() const
Definition: Spectrum.h:219
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
double ComputeMaxScale(MDCMPHelper hSums) const
Represent the ratio between two spectra.
Definition: Ratio.h:8
ReweightableSpectrum fMichelNue
Definition: MichelDecomp.h:192
Spectrum MC_NumuComponent() const override
Base class for the various types of spectrum loader.
const Binning bins
Definition: NumuCC_CPiBin.h:8
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const Cut cut
Definition: exporter_fd.C:30
double ComputeScaleFactor(int EIdx, MDCMPHelper hSums) const
Starting code to run through decomp algebra /////////////////////////////////////////////////////////...
TDirectory * dir
Definition: macro.C:5
REGISTER_LOADFROM("BENDecomp", IDecomp, BENDecomp)
void SaveTo(TDirectory *dir, const std::string &name) const
virtual Spectrum NCTotalComponent() const
Definition: IDecomp.h:18
std::vector< Loaders * > loaders
Definition: syst_header.h:386
double GetSum(ReweightableSpectrum spect, int EIdx) const
Sums given spectrum over NME bins, for a given E Bin.
Standard interface to all decomposition techniques.
Definition: IDecomp.h:13
float fNC[pidBins]
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.
Definition: FileReducer.h:10
void SaveCompPlots(TDirectory *dir, std::string comp, Spectrum dcmp, Spectrum mc) const
const HistAxis MEAxis("N_{Michels}", MEBinning, kNMichels)
Definition: MichelDecomp.h:21
Spectrum MC_NCComponent() const override
std::unique_ptr< IDecomp > LoadFrom< IDecomp >(TDirectory *dir, const std::string &label)
Definition: IDecomp.cxx:53
Double_t sum
Definition: plot.C:31
MDCMPFitResults fFitResult
Definition: MichelDecomp.h:199
Float_t e
Definition: plot.C:35
Spectrum MC_AntiNueComponent() const override
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:114
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:222
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.
Definition: UtilsExt.cxx:115
virtual Spectrum NueComponent() const =0
ReweightableSpectrum fMichelNC
Definition: MichelDecomp.h:190
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