BENDecomp.cxx
Go to the documentation of this file.
2 
5 
7 #include "CAFAna/Cuts/BeamCuts.h"
8 #include "CAFAna/Core/HistAxis.h"
9 #include "CAFAna/Core/Loaders.h"
10 #include "CAFAna/Core/Ratio.h"
13 
23 
24 #include "Utilities/rootlogon.C"
25 
26 #include "TDirectory.h"
27 #include "TObjString.h"
28 #include "TH3.h"
29 #include "TH2.h"
30 #include "TGraph.h"
31 #include "TLegend.h"
32 #include "THStack.h"
33 #include "TCanvas.h"
34 #include "TVectorD.h"
35 #include "TLine.h"
36 
37 #include <cassert>
38 
39 namespace ana
40 {
41  REGISTER_LOADFROM("BENDecomp", IDecomp, BENDecomp);
42 
44  SpectrumLoaderBase& loaderData,
45  const HistAxis& axis,
46  const Cut& cut,
47  const SystShifts& shiftMC,
48  const SystShifts& shiftData,
49  const Var& wei,
50  const HistAxis & axisNumuPi,
51  const Cut& cutNumuPi,
52  const double kaonNormalization
53  )
54  : fData (loaderData, axis, cut, shiftData, wei),
55  fNCTot (loaderMC, axis, cut && kIsNC, shiftMC, wei),
56  fNue (loaderMC, axis, cut && kIsBeamNue&&!kIsAntiNu, shiftMC, wei),
57  fAntiNue (loaderMC, axis, cut && kIsBeamNue&& kIsAntiNu, shiftMC, wei),
58  fNumu (loaderMC, axis, cut && kIsNumuCC &&!kIsAntiNu, shiftMC, wei),
59  fAntiNumu(loaderMC, axis, cut && kIsNumuCC && kIsAntiNu, shiftMC, wei),
60  fTotal (loaderMC, axis, cut && kHasNu, shiftMC, wei),
61  fNotNue (loaderMC, axis, cut && kHasNu && !(kIsBeamNue && !kIsAntiNu), shiftMC, wei),
62 
63  fNueCCFromPi (loaderMC, axis, cut && kIsNueCCFromPiDk, shiftMC, wei),
64  fNueCCFromKa (loaderMC, axis, cut && kIsNueCCFromKaDk, shiftMC, wei),
65  fNueCCFromOther (loaderMC, axis, cut && kIsNueCCNotFromPiDk && kIsNueCCNotFromKaDk, shiftMC, wei),
66  fNueCCFromPiPtPz(loaderMC, ktpzAxis, ktpTAxis, axis, cut && kIsNueCCFromPiDk, shiftMC, wei, Spectrum::kSparse),
67 
68  fNumuSelData (loaderData, axisNumuPi, cutNumuPi, shiftData, wei),
69  fNumuSelBkg (loaderMC, axisNumuPi, cutNumuPi && (!kIsNumuCC || kIsAntiNu), shiftMC, wei),
70  fNumuSelCCFromPi (loaderMC, axisNumuPi, cutNumuPi && kIsNumuCCFromPiDk, shiftMC, wei),
71  fNumuSelCCFromKa (loaderMC, axisNumuPi, cutNumuPi && kIsNumuCCFromKaDk, shiftMC, wei),
72  fNumuSelCCFromOther (loaderMC, axisNumuPi, cutNumuPi && kIsNumuCCNotFromPiDk && kIsNumuCCNotFromKaDk,shiftMC, wei),
73  fNumuSelCCFromPiPtPz(loaderMC, ktpzAxis, ktpTAxis, axisNumuPi, cutNumuPi && kIsNumuCCFromPiDk, shiftMC, wei, Spectrum::kSparse),
74 
75  fNumuUncontainData(Spectrum::Uninitialized()),
76  fNumuUncontainBkg(Spectrum::Uninitialized()),
77  fNumuUncontainCCFromPi(Spectrum::Uninitialized()),
78  fNumuUncontainCCFromKa(Spectrum::Uninitialized()),
79  fNumuUncontainCCFromOther(Spectrum::Uninitialized()),
80 
81  fNueCCFromPiEstim(Spectrum::Uninitialized()),
82  fNueCCFromKaEstim(Spectrum::Uninitialized()),
83 
84  fIsDecomposed(false),
85  fNumuPiWeights(0),
86  fIsFixedKaScale(true),
87  fKaonNormalization(kaonNormalization)
88  {}
89 
91  const HistAxis& axis,
92  const Cut& cut,
93  const SystShifts& shiftMC,
94  const SystShifts& shiftData,
95  const Var& wei,
96  const HistAxis & axisNumuPi,
97  const Cut& cutNumuPi,
98  const double kaonNormalization
99  )
100  : BENDecomp(loaders.GetLoader(caf::kNEARDET, Loaders::kMC),
101  loaders.GetLoader(caf::kNEARDET, Loaders::kData),
102  axis, cut, shiftMC, shiftData, wei,
103  axisNumuPi, cutNumuPi, kaonNormalization)
104  {}
105 
107  SpectrumLoaderBase& loaderMC,
108  SpectrumLoaderBase& loaderData,
109  const HistAxis& axis,
110  const Cut& cut,
111  const SystShifts& shiftMC,
112  const SystShifts& shiftData,
113  const Var& wei,
114  const HistAxis& axisNumuPi,
115  const Cut& cutNumuPi,
116  const HistAxis& axisNumuKa,
117  const Cut& cutNumuKa
118  )
119  : fData (loaderData, axis, cut, shiftData, wei),
120  fNCTot (loaderMC, axis, cut && kIsNC, shiftMC, wei),
121  fNue (loaderMC, axis, cut && kIsBeamNue&&!kIsAntiNu, shiftMC, wei),
122  fAntiNue (loaderMC, axis, cut && kIsBeamNue&& kIsAntiNu, shiftMC, wei),
123  fNumu (loaderMC, axis, cut && kIsNumuCC &&!kIsAntiNu, shiftMC, wei),
124  fAntiNumu(loaderMC, axis, cut && kIsNumuCC && kIsAntiNu, shiftMC, wei),
125  fTotal (loaderMC, axis, cut && kHasNu, shiftMC, wei),
126  fNotNue (loaderMC, axis, cut && !(kIsBeamNue && !kIsAntiNu) && kHasNu, shiftMC, wei),
127 
128  fNueCCFromPi (loaderMC, axis, cut && kIsNueCCFromPiDk, shiftMC, wei),
129  fNueCCFromKa (loaderMC, axis, cut && kIsNueCCFromKaDk, shiftMC, wei),
130  fNueCCFromOther (loaderMC, axis, cut && kIsNueCCNotFromPiDk && kIsNueCCNotFromKaDk, shiftMC, wei),
131  fNueCCFromPiPtPz(loaderMC, ktpzAxis, ktpTAxis, axis, cut && kIsNueCCFromPiDk, shiftMC, wei, Spectrum::kSparse),
132 
133  fNumuSelData (loaderData, axisNumuPi, cutNumuPi, shiftData, wei),
134  fNumuSelBkg (loaderMC, axisNumuPi, cutNumuPi && (!kIsNumuCC || kIsAntiNu), shiftMC, wei),
135  fNumuSelCCFromPi (loaderMC, axisNumuPi, cutNumuPi && kIsNumuCCFromPiDk, shiftMC, wei),
136  fNumuSelCCFromKa (loaderMC, axisNumuPi, cutNumuPi && kIsNumuCCFromKaDk, shiftMC, wei),
137  fNumuSelCCFromOther (loaderMC, axisNumuPi, cutNumuPi && kIsNumuCCNotFromPiDk && kIsNumuCCNotFromKaDk,shiftMC, wei),
138  fNumuSelCCFromPiPtPz(loaderMC, ktpzAxis, ktpTAxis, axisNumuPi, cutNumuPi && kIsNumuCCFromPiDk, shiftMC, wei, Spectrum::kSparse),
139 
140  fNumuUncontainData(loaderData, axisNumuKa, cutNumuKa, shiftData, wei),
141  fNumuUncontainBkg (loaderMC, axisNumuKa, cutNumuKa && (!kIsNumuCC || kIsAntiNu), shiftMC, wei),
142  fNumuUncontainCCFromPi(loaderMC, axisNumuKa, cutNumuKa && kIsNumuCCFromPiDk, shiftMC, wei),
143  fNumuUncontainCCFromKa(loaderMC, axisNumuKa, cutNumuKa && kIsNumuCCFromKaDk, shiftMC, wei),
144  fNumuUncontainCCFromOther(loaderMC, axisNumuKa, cutNumuKa && kIsNumuCCNotFromPiDk && kIsNumuCCNotFromKaDk, shiftMC, wei),
145 
146  fNueCCFromPiEstim(Spectrum::Uninitialized()),
147  fNueCCFromKaEstim(Spectrum::Uninitialized()),
148 
149  fIsDecomposed(false),
150  fNumuPiWeights(0),
151  fIsFixedKaScale(false),
153  {}
154 
156  Loaders& loaders,
157  const HistAxis& axis,
158  const Cut& cut,
159  const SystShifts& shiftMC,
160  const SystShifts& shiftData,
161  const Var& wei,
162  const HistAxis& axisNumuPi,
163  const Cut& cutNumuPi,
164  const HistAxis& axisNumuKa,
165  const Cut& cutNumuKa
166  )
167  : BENDecomp(kaopt ,
168  loaders.GetLoader(caf::kNEARDET, Loaders::kMC),
169  loaders.GetLoader(caf::kNEARDET, Loaders::kData),
170  axis, cut, shiftMC, shiftData, wei,
171  axisNumuPi, cutNumuPi, axisNumuKa, cutNumuKa)
172  {}
173 
174 
175 //--------------------- 38 North -------------------
176 
178  if(!fIsDecomposed)
179  Decompose();
180 
181  double thePOT = fNueCCFromPiEstim.POT();
182 
183  Eigen::ArrayXd aData = fData.GetEigen(thePOT);
184  Eigen::ArrayXd aNue1D = (fNueCCFromPiEstim + fNueCCFromKaEstim + fNueCCFromOther).GetEigen(thePOT);
185  for(int i = 0; i < aNue1D.size(); ++i){
186  // Nue estimate always between 0 and data
187  aNue1D[i] = std::min(std::max(0., aNue1D[i]), aData[i]);
188  }
189 
190  return Spectrum(std::move(aNue1D),
192  thePOT, 0);
193  }
194 
195  void BENDecomp::Decompose() const{
198 
201 
202  fIsDecomposed = true;
203  }
204 
206  {
207  delete fNumuPiWeights;
208  }
209 
211  DontAddDirectory guard;
212 
213  auto thePOT = fNueCCFromPiPtPz.POT();
214  auto hNue3D = fNueCCFromPiPtPz.ToTH3(thePOT); //pt pz var
215  int nbinsx= hNue3D->GetNbinsX();
216  int nbinsy= hNue3D->GetNbinsY();
217  int nbinsz= hNue3D->GetNbinsZ();
218 
219  for (int ee = 1; ee < nbinsz+1; ee++){
220  for (int pz = 1; pz < nbinsx+1; pz++){
221  for (int pt = 1; pt < nbinsy+1; pt++){
222  double ww = fNumuPiWeights->GetBinContent(pz,pt);
223  double nnu = hNue3D->GetBinContent(pz,pt,ee);
224  if(nnu>0 && ww>0) hNue3D->SetBinContent(pz,pt,ee, ww*nnu);
225  }
226  }
227  }
228 
229  std::unique_ptr<TH1D> res(hNue3D->ProjectionZ(UniqueName().c_str()));
230 
231  fNueCCFromPiEstim = Spectrum(Eigen::ArrayXd(Eigen::Map<Eigen::ArrayXd>(res->GetArray(), res->GetNbinsX()+2)),
233  thePOT, 0);
234  delete hNue3D;
235  }
236 
240  }
241 
243  if(!fIsFixedKaScale){
244  double pot = fNumuSelData.POT();
245  double pionScale = 1.; double kaonScale = 1.;
246  double epsilon = 1;
247  //hack hack hack hack
248  //pion contained .75 to 3 GeV
249  //kaon uncontained 4.5 to 10
250 
251  std::unique_ptr<TH1D> hpionContained(fNumuSelCCFromPi.ToTH1(pot));
252  int bin_lo_pi = hpionContained->FindBin(0.75);
253  int bin_hi_pi = hpionContained->FindBin(3);
254  std::unique_ptr<TH1D> hkaonUncontain(fNumuUncontainCCFromKa.ToTH1(pot));
255  int bin_lo_ka = hkaonUncontain->FindBin(4.5);
256  int bin_hi_ka = hkaonUncontain->FindBin(10);
257 
258  while(epsilon>1E-6) {
259 
260  Spectrum kaonCorrect = fNumuSelCCFromKa;
261  kaonCorrect.Scale (kaonScale);
262  Spectrum rawContained = (fNumuSelData - kaonCorrect - fNumuSelBkg -
264 
265  std::unique_ptr<TH1D> hrawContained(rawContained.ToTH1(pot));
266 
267  pionScale = (hrawContained->Integral(bin_lo_pi,bin_hi_pi))/
268  (hpionContained->Integral(bin_lo_pi,bin_hi_pi));
269 
270  Spectrum pionCorrect = fNumuUncontainCCFromPi;
271  pionCorrect.Scale(pionScale);
272 
273  Spectrum rawUncontain = (fNumuUncontainData - pionCorrect - fNumuUncontainBkg -
275 
276 
277  std::unique_ptr<TH1D> hrawUncontain(rawUncontain.ToTH1(pot));
278  epsilon = kaonScale;
279  kaonScale = (hrawUncontain->Integral(bin_lo_ka,bin_hi_ka))/(hkaonUncontain->Integral(bin_lo_ka,bin_hi_ka));
280  epsilon = abs(epsilon - kaonScale)/kaonScale;
281  }
282  fKaonNormalization = kaonScale;
283  }
284 
285 }
286 
288  double thePOT = fNumuSelData.POT();
289  std::unique_ptr<TH1D> hratioNumu(((fNumuSelData - fNumuSelBkg - fNumuSelCCFromOther)/(fNumuSelCCFromPi + fNumuSelCCFromKa)).ToTH1());
290  auto hNumu3D = fNumuSelCCFromPiPtPz.ToTH3(thePOT);
291  TH2* tmp = (TH2*)hNumu3D->Project3D("yx");
292  auto hWeight2D = (TH2*) tmp->Clone();
293  delete tmp;
294  tmp = NULL;
295 
296  int nbinsx= hNumu3D->GetNbinsX();
297  int nbinsy= hNumu3D->GetNbinsY();
298  int nbinsz= hNumu3D->GetNbinsZ();
299 
300  for (int pz = 1; pz < nbinsx+1; pz++){
301  for (int pt = 1; pt < nbinsy+1; pt++){
302  double sum=0,wsum=0/*,tmpn=0*/;
303  for (int ee = 1; ee < nbinsz+1; ee++){
304  double nnu = hNumu3D->GetBinContent(pz,pt,ee);
305  double wei = hratioNumu->GetBinContent(ee);
306  double ene = hratioNumu->GetBinCenter (ee);
307  wsum+= nnu * wei*ene;
308  sum+= nnu * ene;
309 // tmpn+= nnu;
310  }
311  if(wsum > 0) hWeight2D->SetBinContent(pz,pt,wsum/sum);
312 /*
313  if(wsum <= 0 && tmpn > 0) {
314  std::cout<< "WARNING: Check out MakeWeightsNumuFromPion() nnu = "
315  << tmpn <<" while weight not updated here: " <<hWeight2D->GetBinContent(pz, pt); // seems to be ok so far...
316  }
317 */
318  }
319  }
320 
321  Simulation();
322  CornerLabel("#nu-beam");
323  //hWeight2D->SetTitle("Pi weights;pz (GeV/c);pT (GeV/c);weight");
324  fNumuPiWeights = (TH2*) hWeight2D->Clone();
325  delete hWeight2D;
326  delete hNumu3D;
327  }
328  //----------------------------------------------------------------------
330  {
331  if(!fIsDecomposed)
332  Decompose();
333  return (fData - NueEstimate())*(fNumu/fNotNue);
334  }
335 
336  //----------------------------------------------------------------------
338  {
339  if(!fIsDecomposed)
340  Decompose();
341  return (fData - NueEstimate())*(fAntiNumu/fNotNue);
342  }
343 
344  //------------------------- WARNING ------------------------------
345  //----------------------------------------------------------------------
347  {
348  if(!fIsDecomposed)
349  Decompose();
350  return (fData - NueEstimate())*(fNCTot/fNotNue);
351  }
352 
353  //----------------------------------------------------------------------
355  {
356  std::cout << "BENDecomp::NCComponent() not implemented" << std::endl;
357  abort();
358  }
359 
360  //----------------------------------------------------------------------
362  {
363  std::cout << "BENDecomp::NCAntiComponent() not implemented" << std::endl;
364  abort();
365  }
366 
367  //----------------------------------------------------------------------
368 
369  //----------------------------------------------------------------------
371  {
372  if(!fIsDecomposed)
373  Decompose();
374  return NueEstimate();
375  }
376  //----------------------------------------------------------------------
378  {
379  if(!fIsDecomposed)
380  Decompose();
381  return (fData - NueEstimate())*(fAntiNue/fNotNue);
382  }
383 
384  //----------------------------------------------------------------------
385 
386 
387  void BENDecomp::SaveTo(TDirectory* dir, const std::string& name) const
388  {
389  TDirectory* tmp = gDirectory;
390 
391  dir = dir->mkdir(name.c_str()); // switch to subdir
392  dir->cd();
393 
394  TObjString("BENDecomp").Write("type");
395 
396  if(fIsFixedKaScale){
397  TVectorD v(1);v[0]=fKaonNormalization;
398  v.Write("KaonNormalization");
399  }
400  else{
401  fNumuUncontainData.SaveTo(dir, "numu_uncontain_data");
402  fNumuUncontainBkg.SaveTo(dir, "numu_uncontain_bkg");
403  fNumuUncontainCCFromPi.SaveTo(dir, "numuCC_uncontain_from_pi");
404  fNumuUncontainCCFromKa.SaveTo(dir, "numuCC_uncontain_from_ka");
405  fNumuUncontainCCFromOther.SaveTo(dir, "numuCC_uncontain_from_other");
406  }
407  fData.SaveTo(dir, "data_comp");
408  fNCTot.SaveTo(dir, "nc_tot_comp");
409  fNue.SaveTo(dir, "nue_comp");
410  fAntiNue.SaveTo(dir, "antinue_comp");
411  fNumu.SaveTo(dir, "numu_comp");
412  fAntiNumu.SaveTo(dir, "antinumu_comp");
413  fTotal.SaveTo(dir, "total_comp");
414  fNotNue.SaveTo(dir, "not_nue_comp");
415 
416  fNueCCFromPi.SaveTo(dir, "nueCC_from_pi");
417  fNueCCFromKa.SaveTo(dir, "nueCC_from_ka");
418  fNueCCFromOther.SaveTo(dir, "nueCC_from_other");
419  fNueCCFromPiPtPz.SaveTo(dir, "nueCC_from_pi_pt_pz");
420 
421  fNumuSelData.SaveTo(dir, "numu_sel_data");
422  fNumuSelBkg.SaveTo(dir, "numu_sel_bkg");
423  fNumuSelCCFromPi.SaveTo(dir, "numuCC_from_pi");
424  fNumuSelCCFromKa.SaveTo(dir, "numuCC_from_ka");
425  fNumuSelCCFromOther.SaveTo(dir, "numuCC_from_other");
426  fNumuSelCCFromPiPtPz.SaveTo(dir, "numuCC_from_pi_pt_pz");
427 
428  dir->Write();
429  delete dir;
430 
431  tmp->cd();
432  }
433 
434  //----------------------------------------------------------------------
435  std::unique_ptr<BENDecomp>
437  {
438  dir = dir->GetDirectory(name.c_str()); // switch to subdir
439  assert(dir);
440 
441  std::unique_ptr<BENDecomp> ret(new BENDecomp);
442 
443  ret->fData = *Spectrum::LoadFrom(dir, "data_comp");
444  ret->fNCTot = *Spectrum::LoadFrom(dir, "nc_tot_comp");
445  ret->fNue = *Spectrum::LoadFrom(dir, "nue_comp");
446  ret->fAntiNue = *Spectrum::LoadFrom(dir, "antinue_comp");
447  ret->fNumu = *Spectrum::LoadFrom(dir, "numu_comp");
448  ret->fAntiNumu = *Spectrum::LoadFrom(dir, "antinumu_comp");
449  ret->fTotal = *Spectrum::LoadFrom(dir, "total_comp");
450  ret->fNotNue = *Spectrum::LoadFrom(dir, "not_nue_comp");
451 
452  ret->fNueCCFromPi = *Spectrum::LoadFrom(dir, "nueCC_from_pi");
453  ret->fNueCCFromKa = *Spectrum::LoadFrom(dir, "nueCC_from_ka");
454  ret->fNueCCFromOther = *Spectrum::LoadFrom(dir, "nueCC_from_other");
455  ret->fNueCCFromPiPtPz = *Spectrum::LoadFrom(dir, "nueCC_from_pi_pt_pz");
456 
457  ret->fNumuSelData = *Spectrum::LoadFrom(dir, "numu_sel_data");
458  ret->fNumuSelBkg = *Spectrum::LoadFrom(dir, "numu_sel_bkg");
459  ret->fNumuSelCCFromPi = *Spectrum::LoadFrom(dir, "numuCC_from_pi");
460  ret->fNumuSelCCFromKa = *Spectrum::LoadFrom(dir, "numuCC_from_ka");
461  ret->fNumuSelCCFromOther = *Spectrum::LoadFrom(dir, "numuCC_from_other");
462  ret->fNumuSelCCFromPiPtPz = *Spectrum::LoadFrom(dir, "numuCC_from_pi_pt_pz");
463 
464  if(dir->GetListOfKeys()->Contains("KaonNormalization")){
465  auto v = *(TVectorD*)dir->Get("KaonNormalization");
466  ret->fKaonNormalization = v[0];
467  ret->fIsFixedKaScale=true;
468  }
469  else{
470  ret->fNumuUncontainData = *Spectrum::LoadFrom(dir, "numu_uncontain_data");
471  ret->fNumuUncontainBkg = *Spectrum::LoadFrom(dir, "numu_uncontain_bkg");
472  ret->fNumuUncontainCCFromPi = *Spectrum::LoadFrom(dir, "numuCC_uncontain_from_pi");
473  ret->fNumuUncontainCCFromKa = *Spectrum::LoadFrom(dir, "numuCC_uncontain_from_ka");
474  ret->fNumuUncontainCCFromOther = *Spectrum::LoadFrom(dir, "numuCC_uncontain_from_other");
475  }
476 
477  ret->fNumuPiWeights = 0;
478 
479  delete dir;
480 
481  return ret;
482  }
483  void BENDecomp::SavePlotsKa(TDirectory * dir){
484  if(fIsFixedKaScale) {
485  std::cerr << "Using a fixed kaon scale " << fKaonNormalization
486  << ". No plots will be saved " << std::endl;
487  return;
488  }
489  if(!fIsDecomposed) Decompose();
490  rootlogon();
491  auto c1 = new TCanvas(UniqueName().c_str(),"c1",800,500);
492  dir->cd();
493  auto pot = fNumuUncontainData.POT();
494  auto hdata = fNumuUncontainData.ToTH1(pot);
495  hdata->SetMarkerStyle(kFullCircle);
496  hdata->Scale(.001);
497  std::vector<std::pair<Spectrum, int> > sp_col = {{fNumuUncontainBkg,kGray+1},
498  {fNumuUncontainCCFromOther,kGreen+2},
499  {fNumuUncontainCCFromKa,kMagenta -9},
500  {fNumuUncontainCCFromPi,kAzure +4}};
501  auto hs = new THStack();
502  for (auto sp:sp_col){
503  auto hist = sp.first.ToTH1(pot);
504  hist->Scale(.001);
505  hist->SetFillColorAlpha(sp.second,0.3);
506  hs->Add(hist, "hist");
507  }
508 
509  hdata->Draw("e1");
510 
511  hdata->GetYaxis()->SetTitle("10^{3} Events/ 11.0 #times 10^{20} POT");
512  hdata->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
513  hs->Draw("same");
514  hdata->Draw("e1 same");
515 
516  TLegend *leg = new TLegend(0.6, 0.55, 0.89,0.89);
517  leg->SetHeader("Uncontained");
518  TGraph * dummy = new TGraph;
519  dummy->SetLineWidth(3);
520  dummy->SetMarkerStyle(kFullCircle); leg->AddEntry(dummy->Clone(), "#nu_{#mu} data","lpe");
521  dummy->SetMarkerStyle(0);
522  dummy->SetLineColor(kBlack); dummy->SetFillColor(kAzure+4); leg->AddEntry(dummy->Clone(),"#nu_{#mu} from #pi^{#pm}","f");
523  dummy->SetLineColor(kBlack); dummy->SetFillColor(kMagenta-9); leg->AddEntry(dummy->Clone(),"#nu_{#mu} from K^{#pm}","f");
524  dummy->SetLineColor(kBlack); dummy->SetFillColor(kGreen +2 ); leg->AddEntry(dummy->Clone(),"#nu_{#mu} from Other","f");
525  dummy->SetLineColor(kBlack); dummy->SetFillColor(kGray +1 ); leg->AddEntry(dummy->Clone(),"Background","f");
526  leg->Draw();
527  CenterTitles(hdata);
528  Preliminary();
529  CornerLabel("#nu-beam");
530  c1->Write("numuUncontained");
531 
532  auto hratioNueKa = (fNueCCFromKaEstim/fNueCCFromKa).ToTH1();
533  hratioNueKa->GetYaxis()->SetTitle("#nu_{e} Ka Estimate / Ka MC");
534  hratioNueKa->GetYaxis()->SetRangeUser(0, 0.2 + hratioNueKa->GetMaximum());
535  hratioNueKa->Draw("hist");
536  hratioNueKa->Write("nue1Dka");
537 
538  }
539 
540  void BENDecomp::SavePlotsPi(TDirectory * dir) {
541  if(!fIsDecomposed) Decompose();
542  rootlogon();
543  auto c1 = new TCanvas(UniqueName().c_str(),"c1",800,600);
544  dir->cd();
545  auto pot = fNumuSelData.POT();
546  auto hdata2 = fNumuSelData.ToTH1(pot);
547  hdata2->SetMarkerStyle(kFullCircle);
548  hdata2->Scale(.001);
549  std::vector<std::pair<Spectrum, int> > sp_col2 = {{fNumuSelBkg,kGray +1},
550  {fNumuSelCCFromOther,kGreen +2},
551  {fNumuSelCCFromKa,kMagenta -9 },
552  {fNumuSelCCFromPi,kAzure +4}};
553  auto hs2 = new THStack();
554  for (auto sp:sp_col2){
555  auto hist2 = sp.first.ToTH1(pot);
556  hist2->SetFillColorAlpha(sp.second,0.3);
557  hist2->Scale(.001);
558  hs2->Add(hist2, "hist");
559  }
560 
561  hdata2->Draw("e1");
562  hdata2->GetYaxis()->SetTitle("10^{3} Events/ 11.0 #times 10^{20} POT");
563  hdata2->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
564  hs2->Draw("same");
565  hdata2->Draw("e1 same");
566  gPad->Update();
567 
568 
569  TLegend *leg2 = new TLegend(0.6, 0.55, 0.76,0.89);
570  leg2->SetHeader("Contained");
571  TGraph * dummy = new TGraph;
572  dummy->SetLineWidth(3);
573  dummy->SetMarkerStyle(kFullCircle); leg2->AddEntry(dummy->Clone(), "#nu_{#mu} data","lpe");
574  dummy->SetMarkerStyle(0);
575  dummy->SetLineColor(kBlack); dummy->SetFillColor(kAzure +4 ); leg2->AddEntry(dummy->Clone(),"#nu_{#mu} from #pi^{#pm}","f");
576  dummy->SetLineColor(kBlack); dummy->SetFillColor(kMagenta -9 ); leg2->AddEntry(dummy->Clone(),"#nu_{#mu} from K^{#pm}","f");
577  dummy->SetLineColor(kBlack); dummy->SetFillColor(kGreen +2); leg2->AddEntry(dummy->Clone(),"#nu_{#mu} from Other","f");
578  dummy->SetLineColor(kBlack); dummy->SetFillColor(kGray +1); leg2->AddEntry(dummy->Clone(),"Background","f");
579 
580  leg2->Draw();
581  CenterTitles(hdata2);
582  Preliminary();
583  CornerLabel("#nu-beam");
584  c1->Write("numuContained");
585 
586  // Numu weights1D
587  auto hratioNumu = ((fNumuSelData - fNumuSelBkg - fNumuSelCCFromOther)/(fNumuSelCCFromPi + fNumuSelCCFromKa)).ToTH1();
588  hratioNumu->GetYaxis()->SetTitle("#nu_{#mu} (Data-Other) / (MC-Other)");
589 
590  TLine *line = new TLine(0.0,1.0,27.0,1.0);
591  line->SetLineWidth(2);
592  line->SetLineColor(kRed);
593  hratioNumu->Draw("hist");
594  CenterTitles(hratioNumu);
595  Simulation();
596  CornerLabel("#nu-beam");
597  line->Draw("same");
598  c1->Write("numu1Dwgt");
599  // Numu pt pz wgt 2D
600  fNumuPiWeights->GetXaxis()->SetTitle("Ancestor forward momentum p_{z} (GeV/c)");
601  fNumuPiWeights->GetYaxis()->SetTitle("Ancestor trans. momentum p_{t} (GeV/c)");
602  fNumuPiWeights->SetTitle("Pion Weights");
603  fNumuPiWeights->Draw("colz");
605  Simulation();
606  CornerLabel("#nu-beam");
607  c1->Write("numu2Dwgt");
608  // Numu pt pz count 2D
610  hnumu3d->GetXaxis()->SetTitle("Ancestor forward momentum p_{z} (GeV/c)");
611  hnumu3d->GetYaxis()->SetTitle("Ancestor trans. momentum p_{t} (GeV/c)");
612  hnumu3d->GetZaxis()->SetTitle("Events");
613  hnumu3d->Project3D("yx")->Draw("colz");
614  Simulation();
615  CornerLabel("#nu-beam");
616  gPad->Update();
617  c1->Write("numu2D");
618  delete hnumu3d;
619  // Nue pt pz count 2D
620  auto hnue3d = fNueCCFromPiPtPz.ToTH3(fNueCCFromPiPtPz.POT());
621  hnue3d->GetXaxis()->SetTitle("Ancestor forward momentum p_{z} (GeV/c)");
622  hnue3d->GetYaxis()->SetTitle("Ancestor trans. momentum p_{t} (GeV/c)");
623  hnue3d->GetZaxis()->SetTitle("Events");
624  //hnue3d->SetTitle(";Ancestor forward momentum p_{z} (GeV/c);Ancestor trans. momentum p_{t} (GeV/c); Events");
625  // hnue3d->GetZaxis()->Scale(1/1000);
626  hnue3d->Project3D("yx")->Draw("colz");
627  hnue3d->SetTitle("");
628  Simulation();
629  CornerLabel("#nu-beam");
630  gPad->Update();
631  c1->Write("nue2D");
632  delete hnue3d;
633  // Nue Pi estimate/mc ratio
634  auto hratioNuePi = (fNueCCFromPiEstim/fNueCCFromPi).ToTH1();
635  hratioNuePi->GetYaxis()->SetTitle("#nu_{e} Pi Estimate / Pi MC");
636  hratioNuePi->GetYaxis()->SetRangeUser(0.94, 0.2+hratioNuePi->GetMaximum());
637  hratioNuePi->SetTitle("");
638  hratioNuePi->Draw("hist");
639  line->Draw("same");
640  Simulation();
641  CornerLabel("#nu-beam");
642  CenterTitles(hratioNuePi);
643  c1->Write("nue1Dpi");
644 
645  }
646 
647  void BENDecomp::SavePlots(TDirectory * dir) {
648  if(!fIsDecomposed) Decompose();
649  rootlogon();
650  dir->cd();
651  auto c1 = new TCanvas(UniqueName().c_str(),"c1",800,500);
652  double thePOT=fData.POT();
653 
654  auto hs = new THStack();
655  auto hs2 = new THStack();
656 
657  std::vector<std::pair<Spectrum, int> > sp_col = {
658  {fNueCCFromOther,kGray+1},
659  {fNueCCFromKaEstim, kMagenta+2},
660  {fNueCCFromPiEstim, kAzure +4}
661  };
662 
663  std::vector<std::pair<Spectrum, int> > sp_col2 = {
664  {fNueCCFromOther,kGray+1},
665  {fNueCCFromKa, kMagenta+2},
666  {fNueCCFromPi, kAzure+4}
667  };
668 
669  for (auto sp:sp_col){
670  auto hist = sp.first.ToTH1(thePOT, sp.second);
671  hs->Add(hist, "hist");
672  }
673 
674  for (auto sp:sp_col2){
675  auto hist = sp.first.ToTH1(thePOT, sp.second);
676  hist->SetLineStyle(2);
677  hs2->Add(hist, "hist");
678  }
679 // fData.ToTH1(thePOT)->Draw("ep");
680  NueEstimate() .ToTH1(thePOT,kPink+9,1)->Draw("hist");
681  hs->Draw("same");
682  hs2->Draw("same");
683  fNue .ToTH1(thePOT,kPink+9,2)->Draw("hist same");
684  NueEstimate() .ToTH1(thePOT,kPink+9,1)->Draw("hist same");
685 
686  gStyle->SetTextSize(0.03);
687  TLegend * leg = new TLegend(0.7, 0.6, 0.89,0.89);
688  TGraph * dummy = new TGraph;
689  dummy->SetLineWidth(3);
690  dummy->SetLineColor(kPink+9); leg->AddEntry(dummy->Clone(),"Total Beam #nu_{e}","l");
691  dummy->SetLineColor(kAzure+4); leg->AddEntry(dummy->Clone(),"#nu_{e} from Pi","l");
692  dummy->SetLineColor(kMagenta+2); leg->AddEntry(dummy->Clone(),"#nu_{e} from Ka","l");
693  dummy->SetLineColor(kGray+1); leg->AddEntry(dummy->Clone(),"#nu_{e} from Other","l");
694  dummy->SetLineColor(kBlack ); dummy->SetLineStyle(2);
695  leg->AddEntry(dummy->Clone(),"Uncorrected MC","l");
696  leg->Draw();
697  Simulation();
698  CornerLabel("#nu-beam");
699  c1->Write("BENresult");
700 
701  // Nue estimate/mc ratio
702  auto hratioNue = (NueEstimate()/fNue).ToTH1();
703  auto hratioNuePi = (fNueCCFromPiEstim/fNueCCFromPi).ToTH1();
704  auto hratioNueKa = (fNueCCFromKaEstim/fNueCCFromKa).ToTH1();
705  hratioNue->GetYaxis()->SetTitle("#nu_{e} Estimate / MC");
706  hratioNue->GetYaxis()->SetRangeUser(0, hratioNue->GetMaximum()+0.2);
707  hratioNue->SetLineColor(kPink+9);
708  hratioNuePi->SetLineColor(kAzure+4);
709  hratioNueKa->SetLineColor(kMagenta+2);
710  hratioNue->Draw("hist");
711  hratioNue->SetTitle("");
712  hratioNuePi->Draw("hist same");
713  hratioNueKa->Draw("hist same");
714 // leg->GetListOfPrimitives()->RemoveLast();
715 // leg->GetListOfPrimitives()->RemoveLast();
716 // leg->SetY1(0.15);
717 // leg->SetY2(0.45);
718  leg->Draw();
719  CenterTitles(hratioNue);
720  Simulation();
721  CornerLabel("#nu-beam");
722  c1->Write("BENratios");
723  }
724 
725  //----------------------------------------------------------------------
726 
728  return fNue;
729  }
731  return fAntiNue;
732  }
734  return fNumu;
735  }
737  return fAntiNumu;
738  }
740  return fNCTot;
741  }
743  return fData;
744  }
745 
747  const caf::SRProxy* srProxy)
748  {
749 
750  ClearMultiNuInfo(sr);
752 
753  if(sr->hdr.det == caf::kFARDET) return;
754 
756 
757  // remove these completely
758  sr->parent = caf::SRParentBranch();
760  sr->me = caf::SRMichelE();
761 
762  // slc branch
763  // keeping some basic slice info for various cuts
764  unsigned int tempncontplanes = sr->slc.ncontplanes;
765  unsigned int tempnhit = sr->slc.nhit;
766  unsigned int tempfirstplane = sr->slc.firstplane;
767  unsigned int templastplane = sr->slc.lastplane;
768  unsigned int tempncellsfromedge = sr->slc.ncellsfromedge;
769  float tempcalE = sr->slc.calE;
770  caf::SRVector3D tempmeanpos = sr->slc.meanpos;
771  caf::SRVector3D tempboxmax = sr->slc.boxmax;
772  caf::SRVector3D tempboxmin = sr->slc.boxmin;
773 
774  sr->slc = caf::SRSlice();
775 
776  sr->slc.ncontplanes = tempncontplanes;
777  sr->slc.nhit = tempnhit;
778  sr->slc.firstplane = tempfirstplane;
779  sr->slc.lastplane = templastplane;
780  sr->slc.ncellsfromedge = tempncellsfromedge;
781  sr->slc.calE = tempcalE;
782  sr->slc.meanpos = tempmeanpos;
783  sr->slc.boxmin = tempboxmin;
784  sr->slc.boxmax = tempboxmax;
785 
786  // energy branch
787  // keeping only a few numu energy vars
788  float temprecotrkcchadE = sr->energy.numu.recotrkcchadE;
789  float temptrkccE = sr->energy.numu.trkccE;
790  float temptrknonqeE = sr->energy.numu.trknonqeE;
791  float temptrkqeE = sr->energy.numu.trkqeE;
792  float tempndhadcaltranE = sr->energy.numu.ndhadcaltranE;
793  float tempndhadcalcatE = sr->energy.numu.ndhadcalcatE;
794  float tempE = sr->energy.numu.E;
795  float temphadcalE = sr->energy.numu.hadcalE;
796  float temphadtrkE = sr->energy.numu.hadtrkE;
797  float tempndtrkcaltranE = sr->energy.numu.ndtrkcaltranE;
798  float tempndtrkcalcatE = sr->energy.numu.ndtrkcalcatE;
799  float tempndtrkcalactE = sr->energy.numu.ndtrkcalactE;
800  float tempndtrklencat = sr->energy.numu.ndtrklencat;
801  float tempndtrklenact = sr->energy.numu.ndtrklenact;
802 
803  sr->energy.numu = caf::SRNumuEnergy();
804 
805  sr->energy.numu.recotrkcchadE = temprecotrkcchadE;
806  sr->energy.numu.trkccE = temptrkccE;
807  sr->energy.numu.trknonqeE = temptrknonqeE;
808  sr->energy.numu.trkqeE = temptrkqeE;
809  sr->energy.numu.ndhadcaltranE = tempndhadcaltranE;
810  sr->energy.numu.ndhadcalcatE = tempndhadcalcatE;
811  sr->energy.numu.E = tempE;
812  sr->energy.numu.hadcalE = temphadcalE;
813  sr->energy.numu.hadtrkE = temphadtrkE;
814  sr->energy.numu.ndtrkcaltranE = tempndtrkcaltranE;
815  sr->energy.numu.ndtrkcalcatE = tempndtrkcalcatE;
816  sr->energy.numu.ndtrkcalactE = tempndtrkcalactE;
817  sr->energy.numu.ndtrklencat = tempndtrklencat;
818  sr->energy.numu.ndtrklenact = tempndtrklenact;
819 
820  ResetBPFEnergy(sr);
821 
822  // sel branch
823 
824  sr->sel.bpfid = caf::SRBpfId();
825  sr->sel.xnuepid = caf::SRXnue();
826 
827  ResetLEMInfo(sr);
828  ResetLIDInfo(sr);
829  ResetRVPInfo(sr);
830  ResetCosRejInfo(sr);
831 
832  //
833  // track branch
834  // need number of cosmic tracks, therefore we create a dummy vector
835  sr->trk.cosmic.tracks.clear();
836  sr->trk.cosmic.tracks.resize(sr->trk.cosmic.ntracks);
837 
838  // need start and stop coordinates
839  if(sr->trk.kalman.ntracks > 0){
840  std::vector<caf::SRVector3D> temptrackstart;
841  std::vector<caf::SRVector3D> temptrackstop;
842  for(int i = 0; i < (int)sr->trk.kalman.tracks.size(); i++){
843  temptrackstart.push_back(sr->trk.kalman.tracks[i].start);
844  temptrackstop.push_back(sr->trk.kalman.tracks[i].stop);
845  }
846  float temptrackslen = sr->trk.kalman.tracks[0].len;
847 
848  sr->trk.kalman.tracks.clear();
849  sr->trk.kalman.tracks.resize(sr->trk.kalman.ntracks);
850 
851  for(int i = 0; i < (int)sr->trk.kalman.tracks.size(); i++){
852  sr->trk.kalman.tracks[i].start = temptrackstart[i];
853  sr->trk.kalman.tracks[i].stop = temptrackstop[i];
854  }
855  sr->trk.kalman.tracks[0].len = temptrackslen;
856  }
857 
858  // vertex branch
859  // need vtx coordinates
860  if(sr->vtx.elastic.IsValid && sr->vtx.elastic.fuzzyk.npng > 0){
861  float tempvtxX = sr->vtx.elastic.vtx.X();
862  float tempvtxY = sr->vtx.elastic.vtx.Y();
863  float tempvtxZ = sr->vtx.elastic.vtx.Z();
864  float tempnshwlid = sr->vtx.elastic.fuzzyk.nshwlid;
865  float tempnpng = sr->vtx.elastic.fuzzyk.npng;
866  float tempdedx0 = sr->vtx.elastic.fuzzyk.png[0].shwlid.lid.dedx0;
867  int templongestidx = sr->vtx.elastic.fuzzyk.longestidx;
868  float templen = sr->vtx.elastic.fuzzyk.png[templongestidx].len;
869 
870  std::vector<caf::SRVector3D> tempshwlidstart;
871  std::vector<caf::SRVector3D> tempshwlidstop;
872  std::vector<caf::SRVector3D> temppngstart;
873  std::vector<caf::SRCVNParticleResult> tempcvnpart;
874  std::vector<float> tempshwlidcale;
875 
876  for(int pngidx = 0; pngidx < (int)sr->vtx.elastic.fuzzyk.nshwlid; pngidx++){
877  tempshwlidstart.push_back(sr->vtx.elastic.fuzzyk.png[pngidx].shwlid.start);
878  tempshwlidstop.push_back(sr->vtx.elastic.fuzzyk.png[pngidx].shwlid.stop);
879  }
880  for(int pngidx = 0; pngidx < (int)sr->vtx.elastic.fuzzyk.npng; pngidx++){
881  temppngstart.push_back(sr->vtx.elastic.fuzzyk.png[pngidx].start);
882  tempcvnpart.push_back(sr->vtx.elastic.fuzzyk.png[pngidx].cvnpart);
883  tempshwlidcale.push_back(sr->vtx.elastic.fuzzyk.png[pngidx].shwlid.calE);
884  }
885 
886 
887  sr->vtx.elastic.IsValid = true;
888  sr->vtx.elastic.vtx.x = tempvtxX;
889  sr->vtx.elastic.vtx.y = tempvtxY;
890  sr->vtx.elastic.vtx.z = tempvtxZ;
891 
892  sr->vtx.elastic.fuzzyk.npng = tempnpng;
893  sr->vtx.elastic.fuzzyk.nshwlid = tempnshwlid;
894  sr->vtx.elastic.fuzzyk.png.resize(tempnpng);
895  for(int pngidx = 0; pngidx < (int)sr->vtx.elastic.fuzzyk.nshwlid; pngidx++){
896  sr->vtx.elastic.fuzzyk.png[pngidx].shwlid.start = tempshwlidstart[pngidx];
897  sr->vtx.elastic.fuzzyk.png[pngidx].shwlid.stop = tempshwlidstop[pngidx];
898  }
899  for(int pngidx = 0; pngidx < (int)sr->vtx.elastic.fuzzyk.npng; pngidx++){
900  sr->vtx.elastic.fuzzyk.png[pngidx].start = temppngstart[pngidx];
901  sr->vtx.elastic.fuzzyk.png[pngidx].cvnpart = tempcvnpart[pngidx];
902  sr->vtx.elastic.fuzzyk.png[pngidx].shwlid.calE = tempshwlidcale[pngidx];
903  }
904  sr->vtx.elastic.fuzzyk.png[0].shwlid.lid.dedx0 = tempdedx0;
905  sr->vtx.elastic.fuzzyk.longestidx = templongestidx;
906  sr->vtx.elastic.fuzzyk.png[templongestidx].len = templen;
907 
909  }
910 
911  // truth branch
912  if(sr->hdr.ismc == 1 && sr->mc.nnu > 0) {
913 
914  sr->mc.nu[0].rwgt.ppfx.vuniv.clear();
915  sr->mc.nu[0].rwgt.ppfx.nvuniv = 0;
916 
918  }
919 
920  }
921 
922  }
923 
924 }
Spectrum fNumuUncontainData
Definition: BENDecomp.h:188
Det_t det
Detector, ND = 1, FD = 2, NDOS = 3.
Definition: SRHeader.h:28
void Simulation()
Definition: tools.h:16
T max(const caf::Proxy< T > &a, T b)
const std::vector< Binning > & GetBinnings() const
Definition: Spectrum.h:268
const XML_Char * name
Definition: expat.h:151
Spectrum fNueCCFromPiPtPz
Definition: BENDecomp.h:179
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
void MakeWeightsNumuFromPion() const
Definition: BENDecomp.cxx:287
Spectrum fNueCCFromPi
Definition: BENDecomp.h:176
const Cut kIsNumuCCNotFromKaDk
Definition: BeamNueCuts.h:48
Spectrum MC_NCTotalComponent() const override
Definition: BENDecomp.cxx:739
bool fIsFixedKaScale
Definition: BENDecomp.h:199
Far Detector at Ash River.
Definition: SREnums.h:11
A 3-vector with more efficient storage than TVector3.
Definition: SRVector3D.h:14
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
Event ID training variables.
void SavePlotsPi(TDirectory *dir)
Definition: BENDecomp.cxx:540
Spectrum fNumuUncontainCCFromPi
Definition: BENDecomp.h:190
SRHeader hdr
Header branch: run, subrun, etc.
TH1D * hist2
Definition: plotHisto.C:12
const Cut kIsNueCCFromPiDk
Definition: BeamNueCuts.h:44
float ndtrklenact
Near detector – muon track length in active region [cm].
Definition: SRNumuEnergy.h:37
float ndtrkcaltranE
Near detector – muon calorimetric energy in transition plane [GeV].
Definition: SRNumuEnergy.h:40
void ResetCosRejInfo(caf::StandardRecord *sr)
void ClearSecondaryTrackInfo(caf::StandardRecord *sr)
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:209
unsigned int lastplane
last plane
Definition: SRSlice.h:27
unsigned int firstplane
first plane
Definition: SRSlice.h:26
size_t ntracks
Definition: SRKalman.h:23
Spectrum fAntiNumu
Definition: BENDecomp.h:172
float trknonqeE
Track length non-quasielastic neutrino energy [GeV].
Definition: SRNumuEnergy.h:26
Spectrum MC_AntiNueComponent() const override
Definition: BENDecomp.cxx:730
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
SRVector3D boxmax
Maximum coordinates box containing all the hits [cm].
Definition: SRSlice.h:47
std::vector< SRFuzzyKProng > png
Vector of 3D prong objects.
Definition: SRFuzzyK.h:19
SRMichelE me
Michel electron branch.
Proxy for caf::StandardRecord.
Definition: SRProxy.h:1993
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
TH2 * fNumuPiWeights
Definition: BENDecomp.h:198
float Y() const
Definition: SRVector3D.h:33
float ndtrkcalcatE
Near detector – muon calorimetric energy in muon catcher [GeV].
Definition: SRNumuEnergy.h:41
void CornerLabel(std::string &s)
Definition: numu_tools.h:145
SRTrainingBranch training
Extra training information for prototyping PIDs etc.
void NueEstimateFromPi() const
Definition: BENDecomp.cxx:210
const Color_t kMC
Spectrum fNueCCFromKa
Definition: BENDecomp.h:177
const Eigen::ArrayXd & GetEigen() const
NB these don&#39;t have POT scaling. For expert high performance ops only!
Definition: Spectrum.h:200
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 .
OStream cerr
Definition: OStream.cxx:7
Spectrum fNueCCFromKaEstim
Definition: BENDecomp.h:195
void abs(TH1 *hist)
unsigned int longestidx
index of longest prong
Definition: SRFuzzyK.h:23
bool fIsDecomposed
Definition: BENDecomp.h:197
Float_t tmp
Definition: plot.C:36
const HistAxis ktpzAxis("Ancestor target forward momentum tpz (GeV/c)", Binning::Simple(60, 0, 60), kTruetpz)
Definition: HistAxes.h:15
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1481
static std::unique_ptr< BENDecomp > LoadFrom(TDirectory *dir, const std::string &name)
Definition: BENDecomp.cxx:436
bool ismc
data or MC? True if MC
Definition: SRHeader.h:27
Spectrum NCTotalComponent() const override
Definition: BENDecomp.cxx:346
Spectrum fNumuSelCCFromPi
Definition: BENDecomp.h:183
void Scale(double c)
Multiply this spectrum by a constant c.
Definition: Spectrum.cxx:321
Spectrum NCComponent() const override
Definition: BENDecomp.cxx:354
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
Spectrum NueComponent() const override
Definition: BENDecomp.cxx:370
Spectrum fNumuSelCCFromKa
Definition: BENDecomp.h:184
Store BDT variables for the short-baseline oscillation study.
Definition: SRXnue.h:12
Spectrum fTotal
Definition: BENDecomp.h:173
Breakpoint ID (BpfId) output.
Definition: SRBpfId.h:18
Spectrum fData
Definition: BENDecomp.h:167
float X() const
Definition: SRVector3D.h:32
static void ReduceForBEN2020Decaf(caf::StandardRecord *sr, const caf::SRProxy *srProxy)
Definition: BENDecomp.cxx:746
void ResetRVPInfo(caf::StandardRecord *sr)
Spectrum fNumuSelCCFromOther
Definition: BENDecomp.h:185
float trkccE
Track length cc neutrino energy [GeV].
Definition: SRNumuEnergy.h:27
SRKalman kalman
Tracks produced by KalmanTrack.
Definition: SRTrackBranch.h:24
unsigned int nshwlid
number of shwlid showers - either 0 or number of 3d prongs
Definition: SRFuzzyK.h:24
void rootlogon()
Definition: rootlogon.C:7
SRVector3D boxmin
Minimum coordinates box containing all the hits [cm].
Definition: SRSlice.h:46
Spectrum MC_NumuComponent() const override
Definition: BENDecomp.cxx:733
const Cut kIsNumuCCNotFromPiDk
Definition: BeamNueCuts.h:43
Spectrum AntiNueComponent() const override
Definition: BENDecomp.cxx:377
float ndhadcalcatE
Near detector – hadronic calorimetric energy NOT on the muon track in muon catcher [GeV]...
Definition: SRNumuEnergy.h:44
float trkqeE
Track length quasielastic neutrino energy [GeV].
Definition: SRNumuEnergy.h:25
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:607
float hadtrkE
Hadronic calorimetric energy on the muon track[GeV].
Definition: SRNumuEnergy.h:36
const Cut kNue2020NDDecafCut
Definition: NueCuts2020.h:176
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
Float_t E
Definition: plot.C:20
void ClearHoughVertexInfo(caf::StandardRecord *sr)
void ClearMichelTruthInfo(caf::StandardRecord *sr)
Definition: FileReducer.cxx:46
Spectrum fNumuSelData
Definition: BENDecomp.h:181
float calE
Calorimetric energy of the cluster [GeV].
Definition: SRSlice.h:38
const Cut kIsNueCCFromKaDk
Definition: BeamNueCuts.h:49
virtual Spectrum NueEstimate() const
Definition: BENDecomp.cxx:177
#define pot
Numu energy estimator output.
Definition: SRNumuEnergy.h:17
Near Detector in the NuMI cavern.
virtual ~BENDecomp()
Definition: BENDecomp.cxx:205
void SaveTo(TDirectory *dir, const std::string &name) const override
Definition: BENDecomp.cxx:387
Spectrum NumuComponent() const override
Definition: BENDecomp.cxx:329
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:578
SRVector3D vtx
Vertex position in detector coordinates. [cm].
Definition: SRElastic.h:26
short nnu
Number of neutrinos in nu vector (0 or 1)
Definition: SRTruthBranch.h:37
size_t npng
Definition: SRFuzzyK.h:26
float E
Neutrino energy, set to match trkccE [GeV].
Definition: SRNumuEnergy.h:23
An SRSlice contains overarching information for a slice.
Definition: SRSlice.h:15
Spectrum MC_NueComponent() const override
Definition: BENDecomp.cxx:727
float Z() const
Definition: SRVector3D.h:34
unsigned int nhit
number of hits
Definition: SRSlice.h:22
double POT() const
Definition: Spectrum.h:231
const Cut kIsNumuCCFromPiDk
Definition: BeamNueCuts.h:42
void Preliminary()
void MakeWeightsNumuFromKaon() const
Definition: BENDecomp.cxx:242
OStream cout
Definition: OStream.cxx:6
void Decompose() const
Definition: BENDecomp.cxx:195
void NueEstimateFromKa() const
Definition: BENDecomp.cxx:237
SRBpfId bpfid
Output from the BreakPointFitter PID (BPFPIdMaker) package.
Definition: SRIDBranch.h:39
Base class for the various types of spectrum loader.
unsigned int ncellsfromedge
minimum number of cells to edge of detector
Definition: SRSlice.h:30
const Cut kNumuUncontainNDDecafCut
Definition: BeamNueCuts.h:23
Spectrum fNumuUncontainBkg
Definition: BENDecomp.h:189
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
The StandardRecord is the primary top-level object in the Common Analysis File trees.
const Cut cut
Definition: exporter_fd.C:30
void ResetBPFEnergy(caf::StandardRecord *sr)
float ndhadcaltranE
Near detector – hadronic calorimetric energy NOT on the muon track in transition plane [GeV]...
Definition: SRNumuEnergy.h:43
void SavePlots(TDirectory *dir)
Definition: BENDecomp.cxx:647
SRNumuEnergy numu
Numu energy estimator.
const HistAxis ktpTAxis("Ancestor target transv. momentum tpT (GeV/c)", Binning::Simple(10, 0, 1), kTruetpT)
Axes used in Beam Nue Decomposition.
Definition: HistAxes.h:14
void SavePlotsKa(TDirectory *dir)
Definition: BENDecomp.cxx:483
float ndtrkcalactE
Near detector – muon calorimetric energy in active region [GeV].
Definition: SRNumuEnergy.h:39
TH3 * ToTH3(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 3D to obtain TH3.
Definition: Spectrum.cxx:279
Spectrum fNueCCFromOther
Definition: BENDecomp.h:178
void ResetLEMInfo(caf::StandardRecord *sr)
const Cut kHasNu([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return true;})
Definition: TruthCuts.h:56
Int_t nbinsx
Definition: plot.C:23
TDirectory * dir
Definition: macro.C:5
float ndtrklencat
Near detector – muon track length in muon catcher [cm].
Definition: SRNumuEnergy.h:38
REGISTER_LOADFROM("BENDecomp", IDecomp, BENDecomp)
Spectrum AntiNumuComponent() const override
Definition: BENDecomp.cxx:337
Spectrum NCAntiComponent() const override
Definition: BENDecomp.cxx:361
double epsilon
std::vector< Loaders * > loaders
Definition: syst_header.h:386
SRIDBranch sel
Selector (PID) branch.
SRElastic elastic
Single vertex found by Elastic Arms.
GenericHistAxis< Var > HistAxis
Definition: HistAxis.h:111
Spectrum fNotNue
Definition: BENDecomp.h:174
Spectrum fNumuSelCCFromPiPtPz
Definition: BENDecomp.h:186
const Cut kIsNueCCNotFromKaDk
Definition: BeamNueCuts.h:50
assert(nhit_max >=nhit_nbins)
Spectrum fNumuSelBkg
Definition: BENDecomp.h:182
Spectrum fNueCCFromPiEstim
Definition: BENDecomp.h:194
SRSlice slc
Slice branch: nhit, extents, time, etc.
c1
Definition: demo5.py:24
double fKaonNormalization
Definition: BENDecomp.h:200
SRTruthBranch mc
Truth branch for MC: energy, flavor, etc.
SRFuzzyK fuzzyk
Primary 3D prong object.
Definition: SRElastic.h:30
float recotrkcchadE
Reconstructed hadronic energy for track cc neutrino energy estimator [GeV].
Definition: SRNumuEnergy.h:34
This module creates Common Analysis Files.
Definition: FileReducer.h:10
void ResetLIDInfo(caf::StandardRecord *sr)
SRVector3D meanpos
Mean position of hits in slice, weighted by charge [cm].
Definition: SRSlice.h:48
SRParentBranch parent
True parent branch for matching, e.g. MRCC.
Spectrum fNumu
Definition: BENDecomp.h:171
T min(const caf::Proxy< T > &a, T b)
const std::vector< std::string > & GetLabels() const
Definition: Spectrum.h:267
Double_t sum
Definition: plot.C:31
Prevent histograms being added to the current directory.
Definition: Utilities.h:62
const Cut kIsNumuCCFromKaDk
Definition: BeamNueCuts.h:47
unsigned int ncontplanes
number of continuous planes
Definition: SRSlice.h:25
SRTrackBranch trk
Track branch: nhit, len, etc.
Spectrum fNumuUncontainCCFromOther
Definition: BENDecomp.h:192
SREnergyBranch energy
Energy estimator branch.
SRTrackBase cosmic
Tracks produced by CosmicTrack.
Definition: SRTrackBranch.h:26
Spectrum fAntiNue
Definition: BENDecomp.h:170
Spectrum Data_Component() const override
Definition: BENDecomp.cxx:742
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:30
std::vector< SRNeutrino > nu
implemented as a vector to maintain mc.nu structure, i.e. not a pointer, but with 0 or 1 entries...
Definition: SRTruthBranch.h:25
Spectrum fNCTot
Definition: BENDecomp.h:168
std::vector< SRTrack > tracks
Definition: SRTrackBase.h:15
float hadcalE
Hadronic calorimetric energy NOT on the muon track[GeV].
Definition: SRNumuEnergy.h:35
void ClearMultiNuInfo(caf::StandardRecord *sr)
Definition: FileReducer.cxx:29
const Cut kIsNueCCNotFromPiDk
Definition: BeamNueCuts.h:45
static constexpr Double_t sr
Definition: Munits.h:164
std::vector< SRKalmanTrack > tracks
3D Tracks produced by KalmanTrack
Definition: SRKalman.h:16
Spectrum MC_AntiNumuComponent() const override
Definition: BENDecomp.cxx:736
Spectrum fNumuUncontainCCFromKa
Definition: BENDecomp.h:191
const Cut kNumuContainNDDecafCut
Definition: BeamNueCuts.h:22
SRVertexBranch vtx
Vertex branch: location, time, etc.
SRXnue xnuepid
Output from BDT (XnuePID)
Definition: SRIDBranch.h:55
Spectrum fNue
Definition: BENDecomp.h:169