PredictionNueRebinSA.cxx
Go to the documentation of this file.
2 
5 
6 #include "TDirectory.h"
7 #include "TH2.h"
8 #include "TObjString.h"
9 
10 namespace ana
11 {
12  //----------------------------------------------------------------------
14  :
15  fComponentCC{
16  {
17  {pred->ComponentCC(+12, +12), pred->ComponentCC(-12, -12)},
18  {pred->ComponentCC(+12, +14), pred->ComponentCC(-12, -14)},
19  {pred->ComponentCC(+12, +16), pred->ComponentCC(-12, -16)}
20  },
21  {
22  {pred->ComponentCC(+14, +12), pred->ComponentCC(-14, -12)},
23  {pred->ComponentCC(+14, +14), pred->ComponentCC(-14, -14)},
24  {pred->ComponentCC(+14, +16), pred->ComponentCC(-14, -16)}
25  }
26  },
27  fComponentNCTot(pred->ComponentNCTotal())
28  {
29  for(int from = 0; from < 2; ++from)
30  for(int to = 0; to < 3; ++to)
31  for(int anti = 0; anti < 2; ++anti)
32  fComponentCC[from][to][anti] = Reduce(fComponentCC[from][to][anti]);
33 
35  }
36 
37  //----------------------------------------------------------------------
39  {
40  }
41 
42  //----------------------------------------------------------------------
44  {
45  return PredictComponent(calc,
48  Sign::kBoth);
49  }
50 
51  //----------------------------------------------------------------------
53  Flavors::Flavors_t flav,
55  Sign::Sign_t sign) const
56  {
57  Spectrum ret = fComponentNCTot; // Get binning
58  ret.Clear();
59 
60  if(curr & Current::kCC){
61  if(flav & Flavors::kNuEToNuE && sign & Sign::kNu) ret += fComponentCC[0][0][0].Oscillated(calc, +12, +12);
62  if(flav & Flavors::kNuEToNuE && sign & Sign::kAntiNu) ret += fComponentCC[0][0][1].Oscillated(calc, -12, -12);
63 
64  if(flav & Flavors::kNuEToNuMu && sign & Sign::kNu) ret += fComponentCC[0][1][0].Oscillated(calc, +12, +14);
65  if(flav & Flavors::kNuEToNuMu && sign & Sign::kAntiNu) ret += fComponentCC[0][1][1].Oscillated(calc, -12, -14);
66 
67  if(flav & Flavors::kNuEToNuTau && sign & Sign::kNu) ret += fComponentCC[0][2][0].Oscillated(calc, +12, +16);
68  if(flav & Flavors::kNuEToNuTau && sign & Sign::kAntiNu) ret += fComponentCC[0][2][1].Oscillated(calc, -12, -16);
69 
70  if(flav & Flavors::kNuMuToNuE && sign & Sign::kNu) ret += fComponentCC[1][0][0].Oscillated(calc, +14, +12);
71  if(flav & Flavors::kNuMuToNuE && sign & Sign::kAntiNu) ret += fComponentCC[1][0][1].Oscillated(calc, -14, -12);
72 
73  if(flav & Flavors::kNuMuToNuMu && sign & Sign::kNu) ret += fComponentCC[1][1][0].Oscillated(calc, +14, +14);
74  if(flav & Flavors::kNuMuToNuMu && sign & Sign::kAntiNu) ret += fComponentCC[1][1][1].Oscillated(calc, -14, -14);
75 
76  if(flav & Flavors::kNuMuToNuTau && sign & Sign::kNu) ret += fComponentCC[1][2][0].Oscillated(calc, +14, +16);
77  if(flav & Flavors::kNuMuToNuTau && sign & Sign::kAntiNu) ret += fComponentCC[1][2][1].Oscillated(calc, -14, -16);
78  }
79  if(curr & Current::kNC){
80  assert(flav == Flavors::kAll); // Don't know how to calculate anything else
81  assert(sign == Sign::kBoth); // Why would you want to split NCs out by sign?
82 
83  ret += fComponentNCTot;
84  }
85 
86  return Inflate(ret);
87  }
88 
89  //----------------------------------------------------------------------
90  std::unique_ptr<PredictionNueRebinSA> PredictionNueRebinSA::LoadFrom(TDirectory* dir)
91  {
92  TObjString* tag = (TObjString*)dir->Get("type");
93  assert(tag);
94 
95  if(tag->GetString() == "PredictionNueRebinSA"){
96  std::unique_ptr<PredictionNueRebinSA> ret(new PredictionNueRebinSA);
97 
98  for(int i = 0; i < 2; ++i){
99  for(int j = 0; j < 3; ++j){
100  for(int k = 0; k < 2; ++k){
101  ret->fComponentCC[i][j][k] = *ana::LoadFrom<OscillatableSpectrum>(dir->GetDirectory(TString::Format("cc_%d%d%d", i, j, k))).release();
102  }
103  }
104  }
105 
106  ret->fComponentNCTot = *ana::LoadFrom<Spectrum>(dir->GetDirectory("nc")).release();
107 
108  return ret;
109  }
110  else{
111  return std::unique_ptr<PredictionNueRebinSA>(new PredictionNueRebinSA(ana::LoadFrom<IPrediction>(dir).get()));
112  }
113  }
114 
115  //----------------------------------------------------------------------
116  void PredictionNueRebinSA::SaveTo(TDirectory* dir) const
117  {
118  TDirectory* tmp = gDirectory;
119 
120  dir->cd();
121  TObjString("PredictionNueRebinSA").Write("type");
122 
123  for(int i = 0; i < 2; ++i){
124  for(int j = 0; j < 3; ++j){
125  for(int k = 0; k < 2; ++k){
126  fComponentCC[i][j][k].SaveTo(dir->mkdir(TString::Format("cc_%d%d%d", i, j, k)));
127  }
128  }
129  }
130 
131  fComponentNCTot.SaveTo(dir->mkdir("nc"));
132 
133  tmp->cd();
134  }
135 
136  // inbin -> outbin
137  const std::vector<std::pair<int, int>> kBinMapping = {
138  {3, 1}, // 1-1.5 GeV low PID
139  {4, 2},
140  {5, 3},
141  {6, 4}, // 2.5-3 GeV low PID
142  {13, 5}, // 1-1.5 GeV mid PID
143  {14, 6}, // etc
144  {15, 7},
145  {16, 8},
146  {23, 9},
147  {24, 10},
148  {25, 11},
149  {26, 12}
150  };
151 
152  const std::vector<int> kEmptyBins = {0, 1, 2, 7, 8, 9, 10, 11, 12,
153  17, 18, 19, 20, 21, 22,
154  27, 28, 29, 30, 31};
155 
156  //---------------------------------------------------------------------
158  {
159  TH2D* hin = s.ToTH2(1e20);
160  assert(hin->GetNbinsX() == 30); // nue 2D binning
161 
162  // Only 12 bins are really filled
163  TH2D* hout = HistCache::NewTH2D(hin->GetTitle(),
164  Binning::Simple(12, 0, 12));
165 
166  for(const std::pair<int, int>& m: kBinMapping){
167  for(int i = 0; i < hin->GetNbinsY()+2; ++i){
168  hout->SetBinContent(m.second, i, hin->GetBinContent(m.first, i));
169  }
170  }
171 
172  // Check empty bins are really empty
173  for(int i: kEmptyBins){
174  for(int j = 0; j < hin->GetNbinsX()+2; ++j){
175  assert(hin->GetBinContent(i, j) == 0);
176  }
177  }
178 
179  HistCache::Delete(hin);
180 
181  return OscillatableSpectrum(hout, s.GetLabels(), s.GetBinnings(), 1e20, 0);
182  }
183 
184  //---------------------------------------------------------------------
186  {
187  DontAddDirectory guard;
188 
189  TH1D* hin = s.ToTH1(s.POT());
190  assert(hin->GetNbinsX() == 30); // nue 2D binning
191 
192  // Only 12 bins are really filled
193  TH1D* hout = HistCache::New(hin->GetTitle(),
194  Binning::Simple(12, 0, 12));
195 
196  for(const std::pair<int, int>& m: kBinMapping){
197  hout->SetBinContent(m.second, hin->GetBinContent(m.first));
198  }
199 
200  // Check empty bins are really empty
201  for(int i: kEmptyBins){
202  assert(hin->GetBinContent(i) == 0);
203  }
204 
205  HistCache::Delete(hin);
206 
207  return Spectrum(std::unique_ptr<TH1D>(hout), s.GetLabels(), s.GetBinnings(), s.POT(), 0);
208  }
209 
210  //---------------------------------------------------------------------
212  {
213  DontAddDirectory guard;
214 
215  TH1D* hin = s.ToTH1(s.POT());
216  assert(hin->GetNbinsX() == 12); // reduced binning
217 
218  // What external observers are expecting to see
219  TH1D* hout = HistCache::New(hin->GetTitle(),
220  Binning::Simple(30, 0, 30));
221 
222  for(const std::pair<int, int>& m: kBinMapping){
223  // hout->SetBinContent(m.first, hin->GetBinContent(m.second));
224  // This is on the critical path for fits, let's do it the fast way
225  hout->fArray[m.first] = hin->fArray[m.second];
226  }
227 
228  HistCache::Delete(hin);
229 
230  return Spectrum(std::unique_ptr<TH1D>(hout), s.GetLabels(), s.GetBinnings(), s.POT(), 0);
231  }
232 }
Spectrum Oscillated(osc::IOscCalculator *calc, int from, int to) const
static TH2D * NewTH2D(const std::string &title, const Binning &bins)
Definition: HistCache.cxx:80
osc::OscCalculatorDumb calc
Oscillation analysis framework, runs over CAF files outside of ART.
Antineutrinos-only.
Definition: IPrediction.h:50
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
(&#39;beam &#39;)
Definition: IPrediction.h:15
static std::unique_ptr< PredictionNueRebinSA > LoadFrom(TDirectory *dir)
Can load any Prediction or an actual PredictionNueRebinSA.
Spectrum Inflate(const Spectrum &s) const
OscillatableSpectrum Reduce(const OscillatableSpectrum &s) const
void Clear()
Definition: Spectrum.cxx:878
Float_t tmp
Definition: plot.C:36
OscillatableSpectrum fComponentCC[2][3][2]
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
PredictionNueRebinSA()
For LoadFrom.
const XML_Char * s
Definition: expat.h:262
Charged-current interactions.
Definition: IPrediction.h:39
Interactions of both types.
Definition: IPrediction.h:42
std::vector< Binning > GetBinnings() const
std::vector< std::string > GetLabels() const
static void Delete(TH1D *&h)
Definition: HistCache.cxx:215
virtual Spectrum Predict(osc::IOscCalculator *calc) const override
const double j
Definition: BetheBloch.cxx:29
double POT() const
Definition: Spectrum.h:263
Neutrinos-only.
Definition: IPrediction.h:49
def sign(x)
Definition: canMan.py:204
(&#39; survival&#39;)
Definition: IPrediction.h:19
const std::vector< std::pair< int, int > > kBinMapping
std::unique_ptr< IPrediction > LoadFrom< IPrediction >(TDirectory *dir)
Definition: IPrediction.cxx:38
TDirectory * dir
Definition: macro.C:5
virtual OscillatableSpectrum ComponentCC(int from, int to) const
Definition: IPrediction.h:103
Neutral-current interactions.
Definition: IPrediction.h:40
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
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
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
Prevent histograms being added to the current directory.
Definition: Utilities.h:66
Spectrum with true energy information, allowing it to be oscillated
const std::vector< int > kEmptyBins
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:38
virtual void SaveTo(TDirectory *dir) const override
virtual Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
std::vector< Binning > GetBinnings() const
Definition: Spectrum.h:299
void SaveTo(TDirectory *dir) const
void SaveTo(TDirectory *dir) const
Definition: Spectrum.cxx:1029
TH2D * ToTH2(double pot) const