PredictionExtendToPeripheral.cxx
Go to the documentation of this file.
3 
7 #include "CAFAna/Core/Ratio.h"
8 
9 #include "TDirectory.h"
10 #include "TObjString.h"
11 
12 #include "TH2.h"
13 #include "TVectorD.h"
14 
15 #include <iostream>
16 
17 namespace ana
18 {
19  REGISTER_LOADFROM("PredictionExtendToPeripheral", IPrediction, PredictionExtendToPeripheral);
20 
21  //----------------------------------------------------------------------
22 
23  //--- leave this for ana2017 style, i dont wanna update all the old macros in cafana, lets set default as bin3
25  : PredictionExtendToPeripheral(predCore, predNoExtrap, Binning::Simple(27,0,27), mergePeripheral, nbins)
26  {
27  }
28 
29  //----------------------------------------------------------------------
31  : fPredCore(predCore),
32  fPredNoExtrap(predNoExtrap),
33  fBins(bins),
34  fMergePeripheral(mergePeripheral),
35  fNPIDBin(nbins),
51  fIsMerged(false)
52  {
53 
54  if(!fMergePeripheral) std::cout << "WARNING: Not merging peripheral. Check inputs." << std::endl;
55  }
56 
57  //----------------------------------------------------------------------
59  {
60  // We don't necessarily own the predictions, leave them alone
61  }
62 
63  // Reduce 1st, 2nd, 9th for core, 1st, 2nd, 2nd last and last for peripheral
64  //---------------------------------------------------------------------
66  {
67  std::vector<std::string> labels = s.GetLabels();
68  std::vector<Binning> bins = s.GetBinnings();
69  double pot = s.POT();
70  double livetime = s.Livetime();
71 
72  std::vector<int> Bins;
73 
74  for(int i = 1; i < fNBinsO+1; i++){ Bins.push_back(i); }
75  int iter = fNPIDBin;
76 
77  Bins.erase(Bins.end()-2, Bins.end()); // bye peripheral
78  Bins.erase(Bins.begin()+(fNCore),Bins.begin()+(fNCore+2));
79 
80  while(iter-1){ // peripheral is done
81  for(int j: {9, 2, 1}){
82  Bins.erase(Bins.begin() +((iter-2)* 9+j-1));
83  }
84  --iter;
85  }
86 
87  Eigen::MatrixXd ain = s.GetEigen(s.POT());
88  assert(s.GetBinnings()[0].NBins() == fNBinsO && "Merge Peripheral Before Reduce!");
89  Eigen::MatrixXd aout(kTrueEnergyBins.NBins()+2, ReduceBinning.NBins()+2);
90  aout.setZero();
91 
92  int i = 0;
93  for(int x: Bins){
94  ++i;
95  for(int j = 0; j < ain.size(); ++j){
96  aout(j, i) = ain(j, x);
97  }
98  }
99 
100  return OscillatableSpectrum(std::move(aout),
101  HistAxis(labels, bins),
102  pot, livetime);
103 
104  }
105 
106  // One for NC
107  //----------------------------------------------------------------------
109  {
110  std::vector<std::string> labels = s.GetLabels();
111  std::vector<Binning> bins = s.GetBinnings();
112  double pot = s.POT();
113 
114  std::vector<int> Bins;
115 
116  for(int i = 1; i < fNBinsO+1; i++){ Bins.push_back(i); }
117  int iter = fNPIDBin;
118 
119  Bins.erase(Bins.end()-2, Bins.end()); // bye peripheral
120  Bins.erase(Bins.begin()+(fNCore),Bins.begin()+fNCore+2);
121 
122  while(iter-1){ // peripheral is done
123  for(int j: {9, 2, 1}){
124  Bins.erase(Bins.begin() + ((iter-2)* 9+j-1));
125  }
126  --iter;
127  }
128 
129  Eigen::ArrayXd ain = s.GetEigen(pot);
130  assert( ain.size()-2 == fNBinsO && "Merge Peripheral Before Reduce!");
131  Eigen::ArrayXd aout(ReduceBinning.NBins()+2);
132  aout.setZero();
133 
134  for(int x = 0; x < fReduceBin+2; ++x){
135  aout[x]= ain[Bins[x-1]];
136  }
137 
138  std::vector<Binning> tmp_bin = {ReduceBinning};
139  return Spectrum(std::move(aout),
140  HistAxis(s.GetLabels(), tmp_bin),
141  s.POT(), s.Livetime());
142  }
143 
144  //TODO: check following binning in Fit
145  //----------------------------------------------------------------------
147  {
148  assert(!fIsMerged && "Can't call Reduce() twice");
149 
162 
164 
165  fIsMerged = true; // dont collapse the whole universe...
166  }
167 
168  //----------------------------------------------------------------------
170  {
171  return PredictComponent(calc,
174  Sign::kBoth);
175  }
176 
177  //----------------------------------------------------------------------
179  {
180  return PredictComponent(calc,
183  Sign::kBoth);
184  }
185 
186  //----------------------------------------------------------------------
188  {
189  Eigen::ArrayXd a = s.GetEigen(s.POT());
190  Eigen::ArrayXd aret(fNBinsO+2);
191  aret.setZero();
192  for (int i = 1; i <= fNCore; i++){
193  aret[i] = a[i]; // Copy Core
194  }
195 
196  // Merge peripheral
197  for(int i = 0; i <= 8; ++i) aret[fPeriBin] += a[fNBinsI-i];
198 
199  return Spectrum(std::move(aret),
201  s.POT(),
202  s.Livetime());
203  }
204 
205  //----------------------------------------------------------------------
207  {
208  Eigen::MatrixXd aret(kTrueEnergyBins.NBins()+2, fNBinsO+2);
209  aret.setZero();
210  Eigen::MatrixXd aextrap = s.GetEigen(s.POT());
211 
212  for (int binX = 1; binX <= fNCore; binX++){
213  for (int binY = 1; binY <= kTrueEnergyBins.NBins(); binY++){
214  aret(binY, binX) = aextrap(binY, binX);
215  }
216  }
217  // Now, project all the bins in columns 28+ down to the 28th
218  for (int binX = fNCore+1; binX <= fNBinsI; binX++){
219  for (int binY = 1; binY <= kTrueEnergyBins.NBins(); binY++){
220  aret(binY, fPeriBin) += aextrap(binY, binX);
221  }
222  }
223 
224  return OscillatableSpectrum(std::move(aret),
226  s.POT(),
227  s.Livetime());
228  }
229 
230  // Now comes the bookkeeping to stitch the two Predictions together
231  //----------------------------------------------------------------------
233  const Spectrum& NoExtrap) const
234  {
235  // NoExtrap has the nominal FD prediction in it, while Core has the
236  // our decomposed prediction. We want to copy over the decomp weights,
237  // which will just be the Core / NoExtrap bin content. That covers
238  // extrapolation in energy bins.
239  Eigen::ArrayXd arat = Ratio(Core,NoExtrap).GetEigen();
240 
241  const int FirstPeripheralBin = fNCore+1;
242  const int LastPeripheralBin = fNBinsI;
243  const int EBinsPerSelBin = fEnergyBins;
244 
245  for (int bin = FirstPeripheralBin; bin<= LastPeripheralBin; bin++){
246  // Copy decomp ratio from high CVN bin to peripheral bin
247  arat[bin] = arat[bin-EBinsPerSelBin];
248  }
249  Ratio rat(std::move(arat), Core.GetLabels(), Core.GetBinnings());
250 
251  Spectrum extrap = NoExtrap*rat;
252 
253  // Done if we don't want to collapse energy spectrum in peripheral sample
254  if (!fMergePeripheral) return extrap;
255 
256  // Now, we have to merge all energy bins in the peripheral sample
257  // This is only applicable to the Nue2017 result, so assume 28 bins
258  return MergePeripheral(extrap);
259  }
260 
261  //----------------------------------------------------------------------
265  {
266  // NoExtrap has the nominal FD prediction in it, while Core has the
267  // our decomposed prediction. We want to copy over the decomp weights,
268  // which will just be the Core / NoExtrap bin content
269 
270  // Calling Unoscillated() will project down to the predicted reco spectra
271 
272  Spectrum recoCore = Core.Unoscillated();
273 
274  Spectrum recoNoExtrap = NoExtrap.Unoscillated();
275  Eigen::ArrayXd arat = Ratio(recoCore,recoNoExtrap).GetEigen();
276 
277  const int FirstPeripheralBin = fNCore+1;
278  const int LastPeripheralBin = fNBinsO;
279 
280  const int EBinsPerSelBin = fEnergyBins;
281  for (int bin = FirstPeripheralBin; bin<= LastPeripheralBin; bin++){
282  // Copy decomp ratio from high CVN bin to peripheral bin
283  arat[bin] = arat[bin-EBinsPerSelBin];
284  }
285  Ratio rat(std::move(arat), recoCore.GetLabels(), recoCore.GetBinnings());
286  NoExtrap.ReweightToRecoSpectrum(rat*recoNoExtrap);
287 
288  // We're done if we don't want to merge peripheral sample energy dist
289  if (!fMergePeripheral) return NoExtrap;
290 
291  return MergePeripheralOsc(NoExtrap);
292  }
293 
294 
295  //----------------------------------------------------------------------
296  template <typename T>
297  Spectrum
299  Flavors::Flavors_t flav,
301  Sign::Sign_t sign) const
302  {
303  Spectrum ret = fNCTot;
304 
305  ret.Clear();
306 
307  if(curr & Current::kCC){
308 
309  if(flav & Flavors::kNuEToNuE && sign & Sign::kNu)
310  ret += fNuEToNuE.Oscillated(calc, +12, +12);
311  if(flav & Flavors::kNuEToNuE && sign & Sign::kAntiNu)
312  ret += fAntiNuEToAntiNuE.Oscillated(calc, -12, -12);
313  if(flav & Flavors::kNuEToNuMu && sign & Sign::kNu)
314  ret += fNuEToNuMu.Oscillated(calc, +12, +14);
315  if(flav & Flavors::kNuEToNuMu && sign & Sign::kAntiNu)
316  ret += fAntiNuEToAntiNuMu.Oscillated(calc, -12, -14);
317  if(flav & Flavors::kNuEToNuTau && sign & Sign::kNu)
318  ret += fNuEToNuTau.Oscillated(calc, +12, +16);
319  if(flav & Flavors::kNuEToNuTau && sign & Sign::kAntiNu)
320  ret += fAntiNuEToAntiNuTau.Oscillated(calc, -12, -16);
321 
322  if(flav & Flavors::kNuMuToNuE && sign & Sign::kNu)
323  ret += fNuMuToNuE.Oscillated(calc, +14, +12);
324  if(flav & Flavors::kNuMuToNuE && sign & Sign::kAntiNu)
325  ret += fAntiNuMuToAntiNuE.Oscillated(calc, -14, -12);
326  if(flav & Flavors::kNuMuToNuMu && sign & Sign::kNu)
327  ret += fNuMuToNuMu.Oscillated(calc, +14, +14);
328  if(flav & Flavors::kNuMuToNuMu && sign & Sign::kAntiNu)
329  ret += fAntiNuMuToAntiNuMu.Oscillated(calc, -14, -14);
330  if(flav & Flavors::kNuMuToNuTau && sign & Sign::kNu)
331  ret += fNuMuToNuTau.Oscillated(calc, +14, +16);
332  if(flav & Flavors::kNuMuToNuTau && sign & Sign::kAntiNu)
333  ret += fAntiNuMuToAntiNuTau.Oscillated(calc, -14, -16);
334 
335  }
336 
337  if(curr & Current::kNC){
338  assert(flav == Flavors::kAll); // Don't know how to calculate anything else
339 
340  if(sign & Sign::kNu) ret += this->ComponentNC();
341  if(sign & Sign::kAntiNu) ret += this->ComponentNCAnti();
342 
343  }
344 
345  return ret;
346  }
347 
348  template Spectrum
350  Flavors::Flavors_t flav,
352  Sign::Sign_t sign) const;
353  template Spectrum
355  Flavors::Flavors_t flav,
356  Current::Current_t curr,
357  Sign::Sign_t sign) const;
358 
359  //----------------------------------------------------------------------
360  Spectrum
362  Flavors::Flavors_t flav,
363  Current::Current_t curr,
364  Sign::Sign_t sign) const
365  {
366  return _PredictComponent(calc, flav, curr, sign);
367  }
368 
369  //----------------------------------------------------------------------
370  Spectrum
372  Flavors::Flavors_t flav,
373  Current::Current_t curr,
374  Sign::Sign_t sign) const
375  {
376  return _PredictComponent(calc, flav, curr, sign);
377  }
378 
379 
380  //----------------------------------------------------------------------
382  {
383  // beam nues and numu survival come from extending the effect on the high
384  // CVN bin into the peripheral sample
385  if (from == +12 && to == +12) return fNuEToNuE;
386  if (from == +12 && to == +14) return fNuEToNuMu;
387  if (from == +12 && to == +16) return fNuEToNuTau;
388 
389  if (from == +14 && to == +12) return fNuMuToNuE;
390  if (from == +14 && to == +14) return fNuMuToNuMu;
391  if (from == +14 && to == +16) return fNuMuToNuTau;
392 
393  if (from == -12 && to == -12) return fAntiNuEToAntiNuE;
394  if (from == -12 && to == -14) return fAntiNuEToAntiNuMu;
395  if (from == -12 && to == -16) return fAntiNuEToAntiNuTau;
396 
397  if (from == -14 && to == -12) return fAntiNuMuToAntiNuE;
398  if (from == -14 && to == -14) return fAntiNuMuToAntiNuMu;
399  if (from == -14 && to == -16) return fAntiNuMuToAntiNuTau;
400 
401  assert(0 && "Not reached");
402  }
403 
404  //----------------------------------------------------------------------
405  //nc
407  {
408  return fNCTot;
409  }
411  {
412  return fNC;
413  }
415  {
416  return fNCAnti;
417  }
418  //end nc
419  //----------------------------------------------------------------------
421  {
422  TDirectory* tmp = gDirectory;
423 
424  dir = dir->mkdir(name.c_str()); // switch to subdir
425  dir->cd();
426 
427  TObjString("PredictionExtendToPeripheral").Write("type");
428 
429  fPredCore->SaveTo(dir, "predCore");
430  fPredNoExtrap->SaveTo(dir, "predNoExtrap");
431 
432  fBins.SaveTo(dir, "bins");
433 
434  TVectorD vMergePeripheral(1);
435  vMergePeripheral[0] = fMergePeripheral;
436  vMergePeripheral.Write("mergePeripheral");
437 
438  dir->Write();
439  delete dir;
440 
441  tmp->cd();
442  }
443 
444  //----------------------------------------------------------------------
445  std::unique_ptr<PredictionExtendToPeripheral> PredictionExtendToPeripheral::LoadFrom(TDirectory* dir, const std::string& name)
446  {
447  dir = dir->GetDirectory(name.c_str()); // switch to subdir
448  assert(dir);
449 
450  auto predCore = ana::LoadFrom<IPrediction>(dir, "predCore").release();
451  PredictionNoExtrap* predNoExtrap = ana::LoadFrom<PredictionNoExtrap>(dir, "predNoExtrap").release();
452 
453  Binning bins= *Binning::LoadFrom(dir, "bins");
454 
455  bool mergeperipheral = true; // default
456 
457  TVectorD* vMergePeripheral = (TVectorD*)dir->Get("mergePeripheral");
458  // An attempt to be backwards compatable to already generated files
459  if (vMergePeripheral){
460  assert(vMergePeripheral->GetNrows() == 1);
461  mergeperipheral = bool((*vMergePeripheral)[0]);
462  }
463 
464  delete dir;
465 
466  return std::unique_ptr<PredictionExtendToPeripheral>(new PredictionExtendToPeripheral(predCore, predNoExtrap, bins, mergeperipheral));
467  }
468 
469  //----------------------------------------------------------------------
471  {
472  delete fPredCore;
473  delete fPredNoExtrap;
474  }
475 }
static std::unique_ptr< Binning > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Binning.cxx:315
const std::vector< Binning > & GetBinnings() const
Definition: Spectrum.h:268
const XML_Char * name
Definition: expat.h:151
Spectrum _PredictComponent(osc::_IOscCalc< T > *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const
Spectrum ReduceHelperNC(const Spectrum &s) const
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 Eigen::ArrayXd & GetEigen() const
Definition: Ratio.h:42
Antineutrinos-only.
Definition: IPrediction.h:50
bool mergePeripheral
(&#39; appearance&#39;)
Definition: IPrediction.h:18
Spectrum MergePeripheral(const Spectrum &s) const
(&#39;beam &#39;)
Definition: IPrediction.h:15
const Eigen::ArrayXd & GetEigen() const
NB these don&#39;t have POT scaling. For expert high performance ops only!
Definition: Spectrum.h:200
void Clear()
Definition: Spectrum.cxx:433
Float_t tmp
Definition: plot.C:36
virtual void SaveTo(TDirectory *dir, const std::string &name) const
void ReweightToRecoSpectrum(const Spectrum &target)
Recale bins so that Unweighted will return target.
osc::OscCalcDumb calc
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
const std::vector< std::string > & GetLabels() const
const XML_Char * s
Definition: expat.h:262
const int nbins
Definition: cellShifts.C:15
Charged-current interactions.
Definition: IPrediction.h:39
Interactions of both types.
Definition: IPrediction.h:42
const double a
OscillatableSpectrum ReduceHelper(const OscillatableSpectrum &s) const
#define pot
const double j
Definition: BetheBloch.cxx:29
const Binning kTrueEnergyBins
Default true-energy bin edges.
Definition: Binning.h:77
const std::vector< Binning > & GetBinnings() const
std::unique_ptr< IPrediction > LoadFrom< IPrediction >(TDirectory *dir, const std::string &label)
Definition: IPrediction.cxx:28
OscillatableSpectrum MergePeripheralOsc(const OscillatableSpectrum &s) const
float bin[41]
Definition: plottest35.C:14
std::vector< float > Spectrum
Definition: Constants.h:527
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
double POT() const
Definition: Spectrum.h:231
Represent the ratio between two spectra.
Definition: Ratio.h:8
Neutrinos-only.
Definition: IPrediction.h:49
OStream cout
Definition: OStream.cxx:6
(&#39; survival&#39;)
Definition: IPrediction.h:19
const Binning bins
Definition: NumuCC_CPiBin.h:8
OscillatableSpectrum ComponentCC(int from, int to) const override
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
virtual Spectrum Predict(osc::IOscCalc *calc) const override
void NoExtrap(std::string InputString, bool RunProd4Cuts=true, bool RunAllCuts=false, bool Debug=false)
Definition: NoExtrap.C:36
Spectrum Oscillated(osc::IOscCalc *calc, int from, int to) const
TDirectory * dir
Definition: macro.C:5
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Binning.cxx:278
double livetime
Definition: saveFDMCHists.C:21
REGISTER_LOADFROM("BENDecomp", IDecomp, BENDecomp)
int NBins() const
Definition: Binning.h:29
GenericHistAxis< Var > HistAxis
Definition: HistAxis.h:111
Neutral-current interactions.
Definition: IPrediction.h:40
assert(nhit_max >=nhit_nbins)
Eigen::MatrixXd GetEigen(double pot) const
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
virtual Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
OscillatableSpectrum ExtendRecoWeightOscillatable(const OscillatableSpectrum &Core, OscillatableSpectrum NoExtrap) const
Prediction that just uses FD MC, with no extrapolation.
const std::vector< std::string > & GetLabels() const
Definition: Spectrum.h:267
All neutrinos, any flavor.
Definition: IPrediction.h:26
(&#39; appearance&#39;)
Definition: IPrediction.h:16
Spectrum with true energy information, allowing it to be oscillated
static std::unique_ptr< PredictionExtendToPeripheral > LoadFrom(TDirectory *dir, const std::string &name)
Spectrum ExtendRecoWeight(const Spectrum &Core, const Spectrum &NoExtrap) const
virtual void SaveTo(TDirectory *dir, const std::string &name) 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:234
def sign(x)
Definition: canMan.py:197