PredictionExtendToPeripheral.cxx
Go to the documentation of this file.
3 
6 #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 namespace ana
16 {
17  //--- leave this for ana2017 style, i dont wanna update all the old macros in cafana, lets set default as bin3
19  : PredictionExtendToPeripheral(predCore, predNoExtrap, Binning::Simple(27,0,27), mergePeripheral, nbins)
20  {
21  }
22 
23  //----------------------------------------------------------------------
25  : fPredCore(predCore),
26  fPredNoExtrap(predNoExtrap),
27  fBins(bins),
28  fMergePeripheral(mergePeripheral),
29  fNPIDBin(nbins)
30  {
31 
32  if(!fMergePeripheral) std::cout << "WARNING: Not merging peripheral. Check inputs." << std::endl;
33 
36 
39 
42 
45 
48 
51 
55 
56  fIsMerged = false;
57  }
58 
59  //----------------------------------------------------------------------
61  {
62  delete fNCTot;
63  delete fNC;
64  delete fNCAnti;
65  // We don't necessarily own the predictions, leave them alone
66 
67  delete fNuEToNuE;
68  delete fAntiNuEToAntiNuE;
69  delete fNuEToNuMu;
70  delete fAntiNuEToAntiNuMu;
71  delete fNuEToNuTau;
72  delete fAntiNuEToAntiNuTau;
73 
74  delete fNuMuToNuE;
75  delete fAntiNuMuToAntiNuE;
76  delete fNuMuToNuMu;
77  delete fAntiNuMuToAntiNuMu;
78  delete fNuMuToNuTau;
79  delete fAntiNuMuToAntiNuTau;
80  }
81 
82  // Reduce 1st, 2nd, 9th for core, 1st, 2nd, 2nd last and last for peripheral
83  //---------------------------------------------------------------------
85  {
86  std::vector<std::string> labels = s.GetLabels();
87  std::vector<Binning> bins = s.GetBinnings();
88  double pot = s.POT();
89  double livetime = s.Livetime();
90 
91  std::vector<int> Bins;
92 
93  for(int i = 1; i < fNBinsO+1; i++){ Bins.push_back(i); }
94  int iter = fNPIDBin;
95 
96  Bins.erase(Bins.end()-2, Bins.end()); // bye peripheral
97  Bins.erase(Bins.begin()+(fNCore),Bins.begin()+(fNCore+2));
98 
99  while(iter-1){ // peripheral is done
100  for(int j: {9, 2, 1}){
101  Bins.erase(Bins.begin() +((iter-2)* 9+j-1));
102  }
103  --iter;
104  }
105 
106  TH2D* hin = s.ToTH2(s.POT());
107  assert( hin -> GetNbinsX() == fNBinsO && "Merge Peripheral Before Reduce!");
108  TH2D* hout = HistCache::NewTH2D(hin -> GetTitle(), ReduceBinning);
109 
110  int i = 0;
111  for(int x: Bins){
112  ++i;
113  for(int j = 0; j < hin -> GetNbinsY()+2; ++j){
114  hout -> SetBinContent(i, j, hin -> GetBinContent(x, j));
115  }
116  }
117  HistCache::Delete(hin);
118  return OscillatableSpectrum(hout, labels, bins, pot, livetime);
119 
120  }
121 
122  // One for NC
123  //----------------------------------------------------------------------
125  {
126  std::vector<std::string> labels = s.GetLabels();
127  std::vector<Binning> bins = s.GetBinnings();
128  double pot = s.POT();
129 
130  std::vector<int> Bins;
131 
132  for(int i = 1; i < fNBinsO+1; i++){ Bins.push_back(i); }
133  int iter = fNPIDBin;
134 
135  Bins.erase(Bins.end()-2, Bins.end()); // bye peripheral
136  Bins.erase(Bins.begin()+(fNCore),Bins.begin()+fNCore+2);
137 
138  while(iter-1){ // peripheral is done
139  for(int j: {9, 2, 1}){
140  Bins.erase(Bins.begin() + ((iter-2)* 9+j-1));
141  }
142  --iter;
143  }
144 
145  TH1D* hin = s.ToTH1(pot);
146  assert( hin -> GetNbinsX() == fNBinsO && "Merge Peripheral Before Reduce!");
147  TH1D* hout = HistCache::New(hin -> GetTitle(), ReduceBinning);
148 
149  double* fromArr = hin -> GetArray();
150  for(int x = 0; x < fReduceBin+2; ++x){
151  hout -> SetBinContent(x, fromArr[Bins[x-1]]);
152  }
153 
154  HistCache::Delete(hin);
155  std::vector<Binning> tmp_bin = {ReduceBinning};
156  return Spectrum(std::unique_ptr<TH1D>(hout), s.GetLabels(), tmp_bin, s.POT(), s.Livetime());
157  }
158 
159  //TODO: check following binning in Fit
160  //----------------------------------------------------------------------
162  {
163  assert(!fIsMerged && "Can't call Reduce() twice");
164 
177 
179 
180  fIsMerged = true; // dont collapse the whole universe...
181  }
182 
183  //----------------------------------------------------------------------
184  void PredictionExtendToPeripheral::Scale(double from, double to, double scale)
185  {
186  if (from == +12 && to == +12) fNuEToNuE->Scale(scale);
187  else if (from == +12 && to == +14) fNuEToNuMu->Scale(scale);
188  else if (from == +12 && to == +16) fNuEToNuTau->Scale(scale);
189 
190  else if (from == +14 && to == +12) fNuMuToNuE->Scale(scale);
191  else if (from == +14 && to == +14) fNuMuToNuMu->Scale(scale);
192  else if (from == +14 && to == +16) fNuMuToNuTau->Scale(scale);
193 
194  else if (from == -12 && to == -12) fAntiNuEToAntiNuE->Scale(scale);
195  else if (from == -12 && to == -14) fAntiNuEToAntiNuMu->Scale(scale);
196  else if (from == -12 && to == -16) fAntiNuEToAntiNuTau->Scale(scale);
197 
198  else if (from == -14 && to == -12) fAntiNuMuToAntiNuE->Scale(scale);
199  else if (from == -14 && to == -14) fAntiNuMuToAntiNuMu->Scale(scale);
200  else if (from == -14 && to == -16) fAntiNuMuToAntiNuTau->Scale(scale);
201  else if (from == 0 && to == 0){
202  fNCTot->Scale(scale);
203  fNC->Scale(scale);
204  fNCAnti->Scale(scale);
205  }
206  else { std::cout<<"PredictionExtendToPeripheral: this is a joke, right... 0.0"<<std::endl; abort();}
207  }
208 
209  //----------------------------------------------------------------------
211  {
212  return PredictComponent(calc,
215  Sign::kBoth);
216  }
217 
218  //----------------------------------------------------------------------
220  {
222  std::unique_ptr<TH1D> hret(HistCache::New("", bin));
223 
224  TH1D* h = s.ToTH1(s.POT());
225  for (int i = 1; i <= fNCore; i++){
226  hret->SetBinContent(i,h->GetBinContent(i)); // Copy Core;
227  }
228 
229  hret->SetBinContent(fPeriBin, h->Integral(fNBinsI-8, fNBinsI)); // Merge peripheral
231 
232  return new Spectrum(std::move(hret),
233  s.GetLabels(),
234  {bin},
235  s.POT(),
236  s.Livetime());
237  }
238 
239  //----------------------------------------------------------------------
241  {
243  std::unique_ptr<TH2D> hret(HistCache::NewTH2D("", bins));
244  TH2D *hextrap = s.ToTH2(s.POT());
245 
246  for (int binX = 1; binX <= fNCore; binX++){
247  for (int binY = 1; binY < hret->GetNbinsY() +1; binY++){
248  hret->SetBinContent(binX,binY,hextrap->GetBinContent(binX,binY));
249  }
250  }
251  // Now, project all the bins in columns 28+ down to the 28th
252  for (int binX = fNCore+1; binX <= fNBinsI; binX++){
253  for (int binY = 1; binY < hret->GetNbinsY() +1; binY++){
254 
255  hret->SetBinContent(fPeriBin, binY, hret->GetBinContent(fPeriBin,binY)+
256  hextrap->GetBinContent(binX,binY));
257  }
258  }
259  HistCache::Delete(hextrap);
260 
261  return new OscillatableSpectrum(// std::move(hret), // no constructor transferring ownership exists, sadly
262  hret.get(),
263  s.GetLabels(),
264  {bins},
265  s.POT(),
266  s.Livetime());
267  }
268 
269  // Now comes the bookkeeping to stitch the two Predictions together
270  //----------------------------------------------------------------------
272  const Spectrum& NoExtrap) const
273  {
274  // NoExtrap has the nominal FD prediction in it, while Core has the
275  // our decomposed prediction. We want to copy over the decomp weights,
276  // which will just be the Core / NoExtrap bin content. That covers
277  // extrapolation in energy bins.
278  TH1 *hrat = Ratio(Core,NoExtrap).ToTH1();
279 
280  const int FirstPeripheralBin = fNCore+1;
281  const int LastPeripheralBin = fNBinsI;
282  const int EBinsPerSelBin = fEnergyBins;
283 
284  for (int bin = FirstPeripheralBin; bin<= LastPeripheralBin; bin++){
285  // Copy decomp ratio from high CVN bin to peripheral bin
286  hrat->SetBinContent(bin,hrat->GetBinContent(bin-EBinsPerSelBin));
287  }
288  Ratio rat(hrat,"Ana Bin");
289 
290  Spectrum extrap = NoExtrap*rat;
291 
292  // Done if we don't want to collapse energy spectrum in peripheral sample
293  if (!fMergePeripheral)
294  return new Spectrum(extrap);
295 
296  // Now, we have to merge all energy bins in the peripheral sample
297  // This is only applicable to the Nue2017 result, so assume 28 bins
298  return MergePeripheral(extrap);
299  }
300 
301  //----------------------------------------------------------------------
305  {
306  // NoExtrap has the nominal FD prediction in it, while Core has the
307  // our decomposed prediction. We want to copy over the decomp weights,
308  // which will just be the Core / NoExtrap bin content
309 
310  // Calling Unoscillated() will project down to the predicted reco spectra
311 
312  Spectrum recoCore = Core.Unoscillated();
313 
314  Spectrum recoNoExtrap = NoExtrap.Unoscillated();
315  TH1 *hrat = Ratio(recoCore,recoNoExtrap).ToTH1();
316 
317  const int FirstPeripheralBin = fNCore+1;
318  const int LastPeripheralBin = fNBinsO;
319 
320  const int EBinsPerSelBin = fEnergyBins;
321  for (int bin = FirstPeripheralBin; bin<= LastPeripheralBin; bin++){
322  // Copy decomp ratio from high CVN bin to peripheral bin
323  hrat->SetBinContent(bin,hrat->GetBinContent(bin-EBinsPerSelBin));
324  }
325  Ratio rat(hrat,"Ana Bin");
326  NoExtrap.ReweightToRecoSpectrum(rat*recoNoExtrap);
327 
328  // We're done if we don't want to merge peripheral sample energy dist
329  if (!fMergePeripheral)
330  return new OscillatableSpectrum(NoExtrap);
331 
332  return MergePeripheralOsc(NoExtrap);
333  }
334 
335 
336  //----------------------------------------------------------------------
337  Spectrum
339  Flavors::Flavors_t flav,
341  Sign::Sign_t sign) const
342  {
343  Spectrum ret = *fNCTot;
344 
345  ret.Clear();
346 
347  if(curr & Current::kCC){
348 
349  if(flav & Flavors::kNuEToNuE && sign & Sign::kNu)
350  ret += fNuEToNuE->Oscillated(calc, +12, +12);
351  if(flav & Flavors::kNuEToNuE && sign & Sign::kAntiNu)
352  ret += fAntiNuEToAntiNuE->Oscillated(calc, -12, -12);
353  if(flav & Flavors::kNuEToNuMu && sign & Sign::kNu)
354  ret += fNuEToNuMu->Oscillated(calc, +12, +14);
355  if(flav & Flavors::kNuEToNuMu && sign & Sign::kAntiNu)
356  ret += fAntiNuEToAntiNuMu->Oscillated(calc, -12, -14);
357  if(flav & Flavors::kNuEToNuTau && sign & Sign::kNu)
358  ret += fNuEToNuTau->Oscillated(calc, +12, +16);
359  if(flav & Flavors::kNuEToNuTau && sign & Sign::kAntiNu)
360  ret += fAntiNuEToAntiNuTau->Oscillated(calc, -12, -16);
361 
362  if(flav & Flavors::kNuMuToNuE && sign & Sign::kNu)
363  ret += fNuMuToNuE->Oscillated(calc, +14, +12);
364  if(flav & Flavors::kNuMuToNuE && sign & Sign::kAntiNu)
365  ret += fAntiNuMuToAntiNuE->Oscillated(calc, -14, -12);
366  if(flav & Flavors::kNuMuToNuMu && sign & Sign::kNu)
367  ret += fNuMuToNuMu->Oscillated(calc, +14, +14);
368  if(flav & Flavors::kNuMuToNuMu && sign & Sign::kAntiNu)
369  ret += fAntiNuMuToAntiNuMu->Oscillated(calc, -14, -14);
370  if(flav & Flavors::kNuMuToNuTau && sign & Sign::kNu)
371  ret += fNuMuToNuTau->Oscillated(calc, +14, +16);
372  if(flav & Flavors::kNuMuToNuTau && sign & Sign::kAntiNu)
373  ret += fAntiNuMuToAntiNuTau->Oscillated(calc, -14, -16);
374 
375  }
376 
377  if(curr & Current::kNC){
378  assert(flav == Flavors::kAll); // Don't know how to calculate anything else
379 
380  if(sign & Sign::kNu) ret += this->ComponentNC();
381  if(sign & Sign::kAntiNu) ret += this->ComponentNCAnti();
382 
383  }
384 
385  return ret;
386  }
387 
388  //----------------------------------------------------------------------
390  {
391  // beam nues and numu survival come from extending the effect on the high
392  // CVN bin into the peripheral sample
393  if (from == +12 && to == +12) return *fNuEToNuE;
394  if (from == +12 && to == +14) return *fNuEToNuMu;
395  if (from == +12 && to == +16) return *fNuEToNuTau;
396 
397  if (from == +14 && to == +12) return *fNuMuToNuE;
398  if (from == +14 && to == +14) return *fNuMuToNuMu;
399  if (from == +14 && to == +16) return *fNuMuToNuTau;
400 
401  if (from == -12 && to == -12) return *fAntiNuEToAntiNuE;
402  if (from == -12 && to == -14) return *fAntiNuEToAntiNuMu;
403  if (from == -12 && to == -16) return *fAntiNuEToAntiNuTau;
404 
405  if (from == -14 && to == -12) return *fAntiNuMuToAntiNuE;
406  if (from == -14 && to == -14) return *fAntiNuMuToAntiNuMu;
407  if (from == -14 && to == -16) return *fAntiNuMuToAntiNuTau;
408 
409  assert(0 && "Not reached");
410  }
411 
412  //----------------------------------------------------------------------
413  //nc
415  {
416  return *fNCTot;
417  }
419  {
420  return *fNC;
421  }
423  {
424  return *fNCAnti;
425  }
426  //end nc
427  //----------------------------------------------------------------------
429  {
430  TDirectory* tmp = gDirectory;
431 
432  dir->cd();
433 
434  TObjString("PredictionExtendToPeripheral").Write("type");
435 
436  fPredCore->SaveTo(dir->mkdir("predCore"));
437  fPredNoExtrap->SaveTo(dir->mkdir("predNoExtrap"));
438 
439  fBins.SaveTo(dir->mkdir("bins"));
440 
441  TVectorD vMergePeripheral(1);
442  vMergePeripheral[0] = fMergePeripheral;
443  vMergePeripheral.Write("mergePeripheral");
444 
445  tmp->cd();
446  }
447 
448  //----------------------------------------------------------------------
449  std::unique_ptr<PredictionExtendToPeripheral> PredictionExtendToPeripheral::LoadFrom(TDirectory* dir)
450  {
451  assert(dir->GetDirectory("predCore"));
452  assert(dir->GetDirectory("predNoExtrap"));
453  assert(dir->GetDirectory("bins"));
454 
455  auto predCore = ana::LoadFrom<IPrediction>(dir->GetDirectory("predCore")).release();
456  PredictionNoExtrap* predNoExtrap = ana::LoadFrom<PredictionNoExtrap>(dir->GetDirectory("predNoExtrap")).release();
457 
458  TDirectory* subdir = dir->GetDirectory("bins");
459  Binning bins= *Binning::LoadFrom(subdir);
460 
461  bool mergeperipheral = true; // default
462 
463  TVectorD* vMergePeripheral = (TVectorD*)dir->Get("mergePeripheral");
464  // An attempt to be backwards compatable to already generated files
465  if (vMergePeripheral){
466  assert(vMergePeripheral->GetNrows() == 1);
467  mergeperipheral = bool((*vMergePeripheral)[0]);
468  }
469 
470  return std::unique_ptr<PredictionExtendToPeripheral>(new PredictionExtendToPeripheral(predCore, predNoExtrap, bins, mergeperipheral));
471  }
472 }
OscillatableSpectrum ComponentCC(int from, int to) const override
Spectrum Oscillated(osc::IOscCalculator *calc, int from, int to) const
virtual Spectrum ComponentNC() const
Definition: IPrediction.h:111
static TH2D * NewTH2D(const std::string &title, const Binning &bins)
Definition: HistCache.cxx:80
Spectrum ComponentNC() const override
Spectrum * ExtendRecoWeight(const Spectrum &Core, const Spectrum &NoExtrap) const
Spectrum ReduceHelperNC(const Spectrum &s) const
osc::OscCalculatorDumb calc
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Oscillation analysis framework, runs over CAF files outside of ART.
Spectrum ComponentNCTotal() const override
Antineutrinos-only.
Definition: IPrediction.h:50
TH1D * ToTH1(Color_t col=kBlack, Style_t style=kSolid) const
Definition: Ratio.cxx:122
General interface to oscillation calculators.
Definition: FwdDeclare.h:15
(&#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:553
virtual void SaveTo(TDirectory *dir) const
subdir
Definition: cvnie.py:7
(&#39;beam &#39;)
Definition: IPrediction.h:15
void Clear()
Definition: Spectrum.cxx:878
Float_t tmp
Definition: plot.C:36
OscillatableSpectrum * ExtendRecoWeightOscillatable(const OscillatableSpectrum &Core, OscillatableSpectrum NoExtrap) const
Spectrum * MergePeripheral(const Spectrum &s) const
void ReweightToRecoSpectrum(const Spectrum &target)
Recale bins so that Unweighted will return target.
OscillatableSpectrum * MergePeripheralOsc(const OscillatableSpectrum &s) const
void Scale(double c)
Multiply this spectrum by a constant c.
Definition: Spectrum.cxx:734
void Scale(double from, double to, double scale)
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
const XML_Char * s
Definition: expat.h:262
const int nbins
Definition: cellShifts.C:15
Double_t scale
Definition: plot.C:25
Charged-current interactions.
Definition: IPrediction.h:39
virtual void SaveTo(TDirectory *dir) const override
Interactions of both types.
Definition: IPrediction.h:42
std::vector< Binning > GetBinnings() const
std::vector< std::string > GetLabels() const
static std::unique_ptr< PredictionExtendToPeripheral > LoadFrom(TDirectory *dir)
bool mergePeripheral
static void Delete(TH1D *&h)
Definition: HistCache.cxx:215
virtual Spectrum ComponentNCAnti() const
Definition: IPrediction.h:113
OscillatableSpectrum ReduceHelper(const OscillatableSpectrum &s) const
#define pot
void SaveTo(TDirectory *dir) const
Definition: Binning.cxx:189
virtual Spectrum Predict(osc::IOscCalculator *calc) const override
virtual void SaveTo(TDirectory *dir) const override
const double j
Definition: BetheBloch.cxx:29
float bin[41]
Definition: plottest35.C:14
double POT() const
Definition: Spectrum.h:263
Represent the ratio between two spectra.
Definition: Ratio.h:8
Neutrinos-only.
Definition: IPrediction.h:49
def sign(x)
Definition: canMan.py:204
OStream cout
Definition: OStream.cxx:6
(&#39; survival&#39;)
Definition: IPrediction.h:19
cout<< "--"<< endl;for(Int_t iP=1;iP<=hyz->GetNbinsX();iP++){for(Int_t iC=1;iC<=hyz->GetNbinsY();iC++){if(hyv->GetBinContent(iP, iC)>-999){goal_hyv-> SetBinContent(iP, iC,-(dy[iP-1][iC-1]))
virtual Spectrum ComponentNCTotal() const
Definition: IPrediction.h:109
const Binning bins
Definition: NumuCC_CPiBin.h:8
OscillatableSpectrum ComponentCC(int from, int to) const override
std::unique_ptr< IPrediction > LoadFrom< IPrediction >(TDirectory *dir)
Definition: IPrediction.cxx:38
void NoExtrap(std::string InputString, bool RunProd4Cuts=true, bool RunAllCuts=false, bool Debug=false)
Definition: NoExtrap.C:35
TDirectory * dir
Definition: macro.C:5
virtual Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
virtual OscillatableSpectrum ComponentCC(int from, int to) const
Definition: IPrediction.h:103
Neutral-current interactions.
Definition: IPrediction.h:40
Spectrum ComponentNCAnti() const override
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
static TH1D * New(const std::string &title, const Binning &bins)
Definition: HistCache.cxx:24
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
Prediction that just uses FD MC, with no extrapolation.
std::vector< std::string > GetLabels() const
Definition: Spectrum.h:298
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
void Scale(double scale) const
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:38
static std::unique_ptr< Binning > LoadFrom(TDirectory *dir)
Definition: Binning.cxx:221
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:266
std::vector< Binning > GetBinnings() const
Definition: Spectrum.h:299
h
Definition: demo3.py:41
TH2D * ToTH2(double pot) const