FluxDecomp.cxx
Go to the documentation of this file.
2 
4 
7 #include "CAFAna/Cuts/BeamCuts.h"
8 #include "CAFAna/Core/Loaders.h"
10 #include "CAFAna/Core/Ratio.h"
12 
13 #include "TDirectory.h"
14 #include "TObjString.h"
15 #include "TH3.h"
16 #include "TH2.h"
17 #include "TGraph.h"
18 #include "TLegend.h"
19 #include "THStack.h"
20 #include "TCanvas.h"
21 #include "TVectorD.h"
22 #include "TLine.h"
23 
24 #include <cassert>
25 
26 namespace ana
27 {
28  REGISTER_LOADFROM("FluxDecomp", IDecomp, FluxDecomp);
29 
30  //----------------------------------------------------------------------
31 
32  //RHC + FHC
33  FluxDecomp::FluxDecomp (const bool IsRHC,
34  SpectrumLoaderBase& loaderMC,
35  SpectrumLoaderBase& loaderData,
36  const HistAxis& axis,
37  const Cut& cut,
38  const SystShifts& shiftMC,
39  const SystShifts& shiftData,
40  const Var& wei,
41  const HistAxis& axisNumuPi,
42  const Cut& cutNumuPi,
43  const HistAxis& axisNumuKa,
44  const Cut& cutNumuKa
45  )
46  :
47  fIsAntiNu(IsRHC ? kIsAntiNu:!kIsAntiNu),
48  fData (loaderData, axis, cut, shiftData, wei),
49 
50  fNCTot (loaderMC, axis, cut && kIsNC, shiftMC, wei),
51  fNC (loaderMC, axis, cut && kIsNC &&!kIsAntiNu, shiftMC, wei),
52  fNCAnti (loaderMC, axis, cut && kIsNC && kIsAntiNu, shiftMC, wei),
53 
54  fNue (loaderMC, axis, cut && kIsBeamNue&&!kIsAntiNu, shiftMC, wei),
55  fAntiNue (loaderMC, axis, cut && kIsBeamNue&& kIsAntiNu, shiftMC, wei),
56  fNumu (loaderMC, axis, cut && kIsNumuCC &&!kIsAntiNu, shiftMC, wei),
57  fAntiNumu(loaderMC, axis, cut && kIsNumuCC && kIsAntiNu, shiftMC, wei),
58  fTotal (loaderMC, axis, cut && kHasNu, shiftMC, wei),
59 //fhc/rhc difference starts from here
60  fNotNue (loaderMC, axis, cut && !(kIsBeamNue && fIsAntiNu)&& kHasNu, shiftMC, wei), //for prop scale the rest component
61 // optimize truth cuts here for rhc or test kaon flux
62 // revert everything here to anti-nu
63  fNueCCFromPi (loaderMC, axis, cut && kIsBeamNue && kIsAllFromPions && fIsAntiNu, shiftMC, wei),
64  fNueCCFromKa (loaderMC, axis, cut && kIsBeamNue && !kIsAllFromPions&& fIsAntiNu, shiftMC, wei),
65  fNueCCFromOther (loaderMC, axis, cut && kIsBeamNue && !kIsAllFromPions && kIsAllFromPions&& fIsAntiNu, shiftMC, wei),// basically, it means, empty
66  fNueCCFromPiPtPz(loaderMC, ktpzAxis, ktpTAxis, axis, cut && kIsBeamNue && kIsAllFromPions&& fIsAntiNu, shiftMC, wei, Spectrum::kSparse),
67 
68  fNumuSelData (loaderData, axisNumuPi, cutNumuPi, shiftData, wei),
69  fNumuSelBkg (loaderMC, axisNumuPi, cutNumuPi && (!kIsNumuCC || !fIsAntiNu), shiftMC, wei),
70  fNumuSelCCFromPi (loaderMC, axisNumuPi, cutNumuPi && kIsNumuCC && kIsAllFromPions&& fIsAntiNu, shiftMC, wei),
71  fNumuSelCCFromKa (loaderMC, axisNumuPi, cutNumuPi && kIsNumuCC && !kIsAllFromPions&& fIsAntiNu, shiftMC, wei),
72  fNumuSelCCFromOther (loaderMC, axisNumuPi, cutNumuPi && kIsNumuCC && !kIsAllFromPions && kIsAllFromPions&& fIsAntiNu,shiftMC, wei),
73  fNumuSelCCFromPiPtPz(loaderMC, ktpzAxis, ktpTAxis, axisNumuPi, cutNumuPi && kIsNumuCC && kIsAllFromPions&& fIsAntiNu, shiftMC, wei, Spectrum::kSparse),
74 
75  fNumuUncontainData(loaderData, axisNumuKa, cutNumuKa, shiftData, wei),
76  fNumuUncontainBkg (loaderMC, axisNumuKa, cutNumuKa && (!kIsNumuCC ||!fIsAntiNu), shiftMC, wei),
77  fNumuUncontainCCFromPi(loaderMC, axisNumuKa, cutNumuKa && kIsNumuCC && kIsAllFromPions&& fIsAntiNu, shiftMC, wei),
78  fNumuUncontainCCFromKa(loaderMC, axisNumuKa, cutNumuKa && kIsNumuCC && !kIsAllFromPions&& fIsAntiNu, shiftMC, wei),
79  fNumuUncontainCCFromOther(loaderMC, axisNumuKa, cutNumuKa && kIsNumuCC &&!kIsAllFromPions && kIsAllFromPions&& fIsAntiNu, shiftMC, wei),
80 
81  fNueCCFromPiEstim(Spectrum::Uninitialized()),
82  fNueCCFromKaEstim(Spectrum::Uninitialized()),
83 
84  fIsDecomposed(false),
85  fNumuPiWeights(0),
86  fKaonNormalization(1),
87  fIsRHC(IsRHC)
88  {}
89 
90  FluxDecomp::FluxDecomp (const bool IsRHC,
92  const HistAxis& axis,
93  const Cut& cut,
94  const SystShifts& shiftMC,
95  const SystShifts& shiftData,
96  const Var& wei,
97  const HistAxis& axisNumuPi,
98  const Cut& cutNumuPi,
99  const HistAxis& axisNumuKa,
100  const Cut& cutNumuKa
101  )
102  : FluxDecomp(IsRHC,
103  loaders.GetLoader(caf::kNEARDET, Loaders::kMC),
104  loaders.GetLoader(caf::kNEARDET, Loaders::kData),
105  axis, cut, shiftMC, shiftData, wei,
106  axisNumuPi, cutNumuPi, axisNumuKa, cutNumuKa)
107  {}
108 
109 //--------------------- 38 North -------------------
110 
112  if(!fIsDecomposed)
113  Decompose();
114 
115  auto thePOT = fNueCCFromPiEstim.POT();
116 
117  Eigen::ArrayXd aData = fData.GetEigen(thePOT);
118  Eigen::ArrayXd aNue1D = (fNueCCFromPiEstim + fNueCCFromKaEstim + fNueCCFromOther).GetEigen(thePOT);
119  for(int i = 0; i < aNue1D.size(); ++i){
120  // Nue estimate always between 0 and data
121  aNue1D[i] = std::min(std::max(0., aNue1D[i]), aData[i]);
122  }
123 
124  return Spectrum(std::move(aNue1D),
126  thePOT, 0);
127  }
128 
129  void FluxDecomp::Decompose() const{
131 // MakeWeightsNumuFromKaon();
132 
135  fIsDecomposed = true;
136  }
137 
139  {
140  delete fNumuPiWeights;
141  }
142 
144  DontAddDirectory guard;
145 
146  auto thePOT = fNueCCFromPiPtPz.POT();
147  auto hNue3D = fNueCCFromPiPtPz.ToTH3(thePOT); //pt pz var
148  int nbinsx= hNue3D->GetNbinsX();
149  int nbinsy= hNue3D->GetNbinsY();
150  int nbinsz= hNue3D->GetNbinsZ();
151 
152  for (int ee = 1; ee < nbinsz+1; ee++){
153  for (int pz = 1; pz < nbinsx+1; pz++){
154  for (int pt = 1; pt < nbinsy+1; pt++){
155  double ww = fNumuPiWeights->GetBinContent(pz,pt);
156  double nnu = hNue3D->GetBinContent(pz,pt,ee);
157  if(nnu>0 && ww>0) hNue3D->SetBinContent(pz,pt,ee, ww*nnu);
158  }
159  }
160  }
161 
162  std::unique_ptr<TH1D> res(hNue3D->ProjectionZ(UniqueName().c_str()));
163 
164  fNueCCFromPiEstim = Spectrum(Eigen::ArrayXd(Eigen::Map<Eigen::ArrayXd>(res->GetArray(), res->GetNbinsX()+2)),
166  thePOT, 0);
167  delete hNue3D;
168  }
169 /*
170  void FluxDecomp::NueEstimateFromKa() const {
171  double thePOT = fNueCCFromKa.POT();
172  auto hnueK = fNueCCFromKa.ToTH1(thePOT);
173  hnueK->Scale(fKaonNormalization);
174  fNueCCFromKaEstim = Spectrum(hnueK,
175  fData.GetLabels(),
176  fData.GetBinnings(),
177  thePOT, 0);
178  }
179 */
181  double thePOT = fData.POT();
182  Eigen::ArrayXd aNue1Dka = ((fData - fNueCCFromPiEstim)*(fNueCCFromKa/(fTotal-fNueCCFromPi))).GetEigen(thePOT);
183  for (int ee = 1; ee < aNue1Dka.size(); ee++){
184  aNue1Dka[ee] = std::max(0., aNue1Dka[ee]);
185  }
186  fNueCCFromKaEstim = Spectrum(aNue1Dka,
188  thePOT, 0);
189  }
190 
191 // TODO check boundaries
193 
194  double pot = fNumuSelData.POT();
195  double pionScale = 1.; double kaonScale = 1.;
196  double epsilon = 1;
197  //hack hack hack hack
198  //pion contained .75 to 3 GeV
199  //kaon uncontained 4.5 to 10
200 
201  auto hpionContained = fNumuSelCCFromPi.ToTH1(pot);
202  int bin_lo_pi = hpionContained->FindBin(0.75);
203  int bin_hi_pi = hpionContained->FindBin(3);
204  auto hkaonUncontain = fNumuUncontainCCFromKa.ToTH1(pot);
205  int bin_lo_ka = hkaonUncontain->FindBin(4.5);
206  int bin_hi_ka = hkaonUncontain->FindBin(10);
207 
208  while(epsilon>1E-6) {
209 
210  Spectrum kaonCorrect = fNumuSelCCFromKa;
211  kaonCorrect.Scale (kaonScale);
212  Spectrum rawContained = (fNumuSelData - kaonCorrect - fNumuSelBkg -
214 
215  auto hrawContained = rawContained.ToTH1(pot);
216 
217  pionScale = (hrawContained->Integral(bin_lo_pi,bin_hi_pi))/
218  (hpionContained->Integral(bin_lo_pi,bin_hi_pi));
219 
220  Spectrum pionCorrect = fNumuUncontainCCFromPi;
221  pionCorrect.Scale(pionScale);
222 
223  Spectrum rawUncontain = (fNumuUncontainData - pionCorrect - fNumuUncontainBkg -
225 
226 
227  auto hrawUncontain = rawUncontain.ToTH1(pot);
228  epsilon = kaonScale;
229  kaonScale = (hrawUncontain->Integral(bin_lo_ka,bin_hi_ka))/(hkaonUncontain->Integral(bin_lo_ka,bin_hi_ka));
230  epsilon = abs(epsilon - kaonScale)/kaonScale;
231  }
232  fKaonNormalization = kaonScale;
233 
234  }
235 
237  double thePOT = fNumuSelData.POT();
239  auto hNumu3D = fNumuSelCCFromPiPtPz.ToTH3(thePOT);
240  TH2* tmp = (TH2*)hNumu3D->Project3D("yx");
241  auto hWeight2D = (TH2*) tmp->Clone();
242  delete tmp;
243 
244  int nbinsx= hNumu3D->GetNbinsX();
245  int nbinsy= hNumu3D->GetNbinsY();
246  int nbinsz= hNumu3D->GetNbinsZ();
247 
248  for (int pz = 1; pz < nbinsx+1; pz++){
249  for (int pt = 1; pt < nbinsy+1; pt++){
250  double sum=0,wsum=0;
251  for (int ee = 1; ee < nbinsz+1; ee++){
252  double nnu = hNumu3D->GetBinContent(pz,pt,ee);
253  double wei = hratioNumu->GetBinContent(ee);
254  double ene = hratioNumu->GetBinCenter (ee);
255  wsum+= nnu * wei*ene;
256  sum+= nnu * ene;
257  }
258  if(wsum > 0) hWeight2D->SetBinContent(pz,pt,wsum/sum);
259  }
260  }
261  hWeight2D->SetTitle("#bar{#nu}_{#mu} from Pi weights;pz (GeV/c);pT (GeV/c);weight");
262  fNumuPiWeights = hWeight2D;
263  delete hNumu3D;
264  }
265  //----------------------------------------------------------------------
267  {
268  if(!fIsDecomposed)
269  Decompose();
270  return (fData - NueEstimate())*(fNumu/fNotNue);
271  }
272 
273  //----------------------------------------------------------------------
275  {
276  if(!fIsDecomposed)
277  Decompose();
278  return (fData - NueEstimate())*(fAntiNumu/fNotNue);
279  }
280 
281  //------------------------- WARNING ------------------------------
282  //----------------------------------------------------------------------
284  {
285  if(!fIsDecomposed)
286  Decompose();
287  return (fData - NueEstimate())*(fNCTot/fNotNue);
288  }
289 
290  //----------------------------------------------------------------------
292  {
293  if(!fIsDecomposed)
294  Decompose();
295  return (fData - NueEstimate())*(fNC/fNotNue);
296 
297  }
298 
299  //----------------------------------------------------------------------
301  {
302  if(!fIsDecomposed)
303  Decompose();
304  return (fData - NueEstimate())*(fNCAnti/fNotNue);
305 
306  }
307 
308  //----------------------------------------------------------------------
309 
310  //----------------------------------------------------------------------
312  {
313  if(!fIsDecomposed)
314  Decompose();
315  if(fIsRHC) return (fData - NueEstimate())*(fNue/fNotNue);// Nue is prop scaled in RHC BEN now:)
316  else return NueEstimate();
317  }
318 //AntiNueComponent is different as in bkg or Nue... nu/nubar ppfx beamnue recoE ratio checked for 2018Ana, F/N not yet
319 //
320  //----------------------------------------------------------------------
322  {
323  if(!fIsDecomposed)
324  Decompose();
325  if(fIsRHC) return NueEstimate();
326  else return (fData - NueEstimate())*(fAntiNue/fNotNue);
327  }
328 
329  //----------------------------------------------------------------------
330  void FluxDecomp::SaveTo(TDirectory* dir, const std::string& name) const
331  {
332  TDirectory* tmp = gDirectory;
333 
334  dir = dir->mkdir(name.c_str()); // switch to subdir
335  dir->cd();
336 
337  TObjString("FluxDecomp").Write("type");
338 
339  if(fIsRHC){
340  TVectorD r(1);
341  r[0] = fIsRHC;
342  r.Write("IsRHC");
343  }
344  fNumuUncontainData.SaveTo(dir, "numu_uncontain_data");
345  fNumuUncontainBkg.SaveTo(dir, "numu_uncontain_bkg");
346  fNumuUncontainCCFromPi.SaveTo(dir, "numuCC_uncontain_from_pi");
347  fNumuUncontainCCFromKa.SaveTo(dir, "numuCC_uncontain_from_ka");
348  fNumuUncontainCCFromOther.SaveTo(dir, "numuCC_uncontain_from_other");
349 
350 
351  fData.SaveTo(dir, "data_comp");
352  fNCTot.SaveTo(dir, "nc_tot_comp");
353  fNC.SaveTo(dir, "nc_comp");
354  fNCAnti.SaveTo(dir, "nc_anti_comp");
355 
356  fNue.SaveTo(dir, "nue_comp");
357  fAntiNue.SaveTo(dir, "antinue_comp");
358  fNumu.SaveTo(dir, "numu_comp");
359  fAntiNumu.SaveTo(dir, "antinumu_comp");
360  fTotal.SaveTo(dir, "total_comp");
361  fNotNue.SaveTo(dir, "not_nue_comp");
362 
363  fNueCCFromPi.SaveTo(dir, "nueCC_from_pi");
364  fNueCCFromKa.SaveTo(dir, "nueCC_from_ka");
365  fNueCCFromOther.SaveTo(dir, "nueCC_from_other");
366  fNueCCFromPiPtPz.SaveTo(dir, "nueCC_from_pi_pt_pz");
367 
368  fNumuSelData.SaveTo(dir, "numu_sel_data");
369  fNumuSelBkg.SaveTo(dir, "numu_sel_bkg");
370  fNumuSelCCFromPi.SaveTo(dir, "numuCC_from_pi");
371  fNumuSelCCFromKa.SaveTo(dir, "numuCC_from_ka");
372  fNumuSelCCFromOther.SaveTo(dir, "numuCC_from_other");
373  fNumuSelCCFromPiPtPz.SaveTo(dir, "numuCC_from_pi_pt_pz");
374 
375  dir->Write();
376  delete dir;
377 
378  tmp->cd();
379  }
380 
381  //----------------------------------------------------------------------
382  std::unique_ptr<FluxDecomp>
384  {
385  dir = dir->GetDirectory(name.c_str()); // switch to subdir
386  assert(dir);
387 
388  std::unique_ptr<FluxDecomp> ret(new FluxDecomp);
389 
390  ret->fData = *Spectrum::LoadFrom(dir, "data_comp");
391  ret->fNCTot = *Spectrum::LoadFrom(dir, "nc_tot_comp");
392  ret->fNC = *Spectrum::LoadFrom(dir, "nc_comp");
393  ret->fNCAnti = *Spectrum::LoadFrom(dir, "nc_anti_comp");
394 
395  ret->fNue = *Spectrum::LoadFrom(dir, "nue_comp");
396  ret->fAntiNue = *Spectrum::LoadFrom(dir, "antinue_comp");
397  ret->fNumu = *Spectrum::LoadFrom(dir, "numu_comp");
398  ret->fAntiNumu = *Spectrum::LoadFrom(dir, "antinumu_comp");
399  ret->fTotal = *Spectrum::LoadFrom(dir, "total_comp");
400  ret->fNotNue = *Spectrum::LoadFrom(dir, "not_nue_comp");
401 
402  ret->fNueCCFromPi = *Spectrum::LoadFrom(dir, "nueCC_from_pi");
403  ret->fNueCCFromKa = *Spectrum::LoadFrom(dir, "nueCC_from_ka");
404  ret->fNueCCFromOther = *Spectrum::LoadFrom(dir, "nueCC_from_other");
405  ret->fNueCCFromPiPtPz = *Spectrum::LoadFrom(dir, "nueCC_from_pi_pt_pz");
406 
407  ret->fNumuSelData = *Spectrum::LoadFrom(dir, "numu_sel_data");
408  ret->fNumuSelBkg = *Spectrum::LoadFrom(dir, "numu_sel_bkg");
409  ret->fNumuSelCCFromPi = *Spectrum::LoadFrom(dir, "numuCC_from_pi");
410  ret->fNumuSelCCFromKa = *Spectrum::LoadFrom(dir, "numuCC_from_ka");
411  ret->fNumuSelCCFromOther = *Spectrum::LoadFrom(dir, "numuCC_from_other");
412  ret->fNumuSelCCFromPiPtPz = *Spectrum::LoadFrom(dir, "numuCC_from_pi_pt_pz");
413 
414  if(dir->GetListOfKeys()->Contains("IsRHC")) ret->fIsRHC = true;
415  else ret->fIsRHC = false;
416 
417  ret->fNumuUncontainData = *Spectrum::LoadFrom(dir, "numu_uncontain_data");
418  ret->fNumuUncontainBkg = *Spectrum::LoadFrom(dir, "numu_uncontain_bkg");
419  ret->fNumuUncontainCCFromPi = *Spectrum::LoadFrom(dir, "numuCC_uncontain_from_pi");
420  ret->fNumuUncontainCCFromKa = *Spectrum::LoadFrom(dir, "numuCC_uncontain_from_ka");
421  ret->fNumuUncontainCCFromOther = *Spectrum::LoadFrom(dir, "numuCC_uncontain_from_other");
422 
423  ret->fNumuPiWeights = 0;
424 
425  delete dir;
426 
427  return ret;
428  }
429  void FluxDecomp::SavePlotsKa(TDirectory * dir){
430 
431  if(fIsRHC) {
432  std::cout<< "Hello CAFAna I'm RHC~~~" << std::endl;
433  }
434 
435  if(!fIsDecomposed) Decompose();
436  auto c1 = new TCanvas(UniqueName().c_str(),"c1",800,500);
437  dir->cd();
438  auto pot = fNumuUncontainData.POT();
439  auto hdata = fNumuUncontainData.ToTH1(pot);
440  hdata->SetMarkerStyle(kFullCircle);
441 
442  std::vector<std::pair<Spectrum, int> > sp_col = {{fNumuUncontainBkg,kGray+1},
443  {fNumuUncontainCCFromOther,kGreen+2},
444  {fNumuUncontainCCFromKa,kMagenta -9},
445  {fNumuUncontainCCFromPi,kAzure +6}};
446  auto hs = new THStack();
447  for (auto sp:sp_col){
448  auto hist = sp.first.ToTH1(pot);
449  hist->SetFillColorAlpha(sp.second,0.3);
450  hs->Add(hist, "hist");
451  }
452 
453  hdata->Draw("e1");
454  hs->Draw("same");
455  hdata->Draw("e1 same");
456  TLegend* leg = new TLegend(0.6, 0.6, 0.89,0.89);
457  leg->SetHeader("Uncontained");
458  TGraph * dummy = new TGraph;
459  dummy->SetLineWidth(3);
460  dummy->SetMarkerStyle(kFullCircle); leg->AddEntry(dummy->Clone(), "#bar{#nu}_{#mu} data","lpe");
461  dummy->SetMarkerStyle(0);
462  dummy->SetLineColor(kAzure+6); leg->AddEntry(dummy->Clone(),"#bar{#nu}_{#mu} from Pi","l");
463  dummy->SetLineColor(kMagenta-9); leg->AddEntry(dummy->Clone(),"#bar{#nu}_{#mu} from Ka","l");
464  dummy->SetLineColor(kGreen +2 ); leg->AddEntry(dummy->Clone(),"#bar{#nu}_{#mu} from Other","l");
465  dummy->SetLineColor(kGray +1 ); leg->AddEntry(dummy->Clone(),"Bkg.","l");
466  leg->Draw();
467  c1->Write("numuUncontained");
468 
469  dir->cd();
470  pot = fNumuSelData.POT();
471  hdata = fNumuSelData.ToTH1(pot);
472  hdata->SetMarkerStyle(kFullCircle);
473 
474  sp_col = {{fNumuSelBkg,kGray +1},
475  {fNumuSelCCFromOther,kGreen +2},
476  {fNumuSelCCFromKa,kMagenta -9 },
477  {fNumuSelCCFromPi,kAzure +6}};
478  hs = new THStack();
479  for (auto sp:sp_col){
480  auto hist = sp.first.ToTH1(pot);
481  hist->SetFillColorAlpha(sp.second,0.3);
482  hs->Add(hist, "hist");
483  }
484 
485  hdata->Draw("e1");
486  hs->Draw("same");
487  hdata->Draw("e1 same");
488  leg = new TLegend(0.6, 0.6, 0.89,0.89);
489  leg->SetHeader("Contained");
490  dummy = new TGraph;
491  dummy->SetLineWidth(3);
492  dummy->SetMarkerStyle(kFullCircle); leg->AddEntry(dummy->Clone(), "#bar{#nu}_{#mu} data","lpe");
493  dummy->SetMarkerStyle(0);
494  dummy->SetLineColor(kAzure +6 ); leg->AddEntry(dummy->Clone(),"#bar{#nu}_{#mu} from Pi","l");
495  dummy->SetLineColor(kMagenta -9 ); leg->AddEntry(dummy->Clone(),"#bar{#nu}_{#mu} from Ka","l");
496  dummy->SetLineColor(kGreen +2); leg->AddEntry(dummy->Clone(),"#bar{#nu}_{#mu} from Other","l");
497  dummy->SetLineColor(kGray +1 ); leg->AddEntry(dummy->Clone(),"Bkg.","l");
498  leg->Draw();
499  c1->Write("numuContained");
500 
501  auto hratioNueKa = (fNueCCFromKaEstim/fNueCCFromKa).ToTH1();
502  hratioNueKa->GetYaxis()->SetTitle("#bar{#nu}_{e} Ka Estimate / Ka MC");
503  TLine *line = new TLine(0.0,1.0,27.0,1.0);
504  line->SetLineWidth(2);
505  line->SetLineColor(kRed);
506  hratioNueKa->Draw("hist");
507  line->Draw("same");
508  c1->Write("nue1Dka");
509 
510  }
511 
512  void FluxDecomp::SavePlotsPi(TDirectory * dir) {
513  if(!fIsDecomposed) Decompose();
514  auto c1 = new TCanvas(UniqueName().c_str(),"c1",800,500);
515  dir->cd();
516  // Numu weights1D
517  auto hratioNumu = ((fNumuSelData - fNumuSelBkg - fNumuSelCCFromOther)/(fNumuSelCCFromPi + fNumuSelCCFromKa)).ToTH1();
518  hratioNumu->GetYaxis()->SetTitle("#bar{#nu}_{#mu} (Data-Other) / (MC-Other)");
519  TLine *line = new TLine(0.0,1.0,27.0,1.0);
520  line->SetLineWidth(2);
521  line->SetLineColor(kRed);
522  hratioNumu->Draw("hist");
523  line->Draw("same");
524  c1->Write("numu1Dwgt");
525  // Numu pt pz wgt 2D
526  fNumuPiWeights->Draw("colz");
527  c1->Write("numu2Dwgt");
528  // Numu pt pz count 2D
530  hnumu3d->SetTitle("#bar{#nu}_{#mu} from Pi;pz (GeV/c);pT (GeV/c); Events");
531  hnumu3d->Project3D("yx")->Draw("colz");
532  c1->Write("numu2D");
533  // Nue pt pz count 2D
534  auto hnue3d = fNueCCFromPiPtPz.ToTH3(fNueCCFromPiPtPz.POT());
535  hnue3d->SetTitle("#bar{#nu}_{e} from Pi;pz (GeV/c);pT (GeV/c); Events");
536  hnue3d->Project3D("yx")->Draw("colz");
537  c1->Write("nue2D");
538  // Nue Pi estimate/mc ratio
539  auto hratioNuePi = (fNueCCFromPiEstim/fNueCCFromPi).ToTH1();
540  hratioNuePi->GetYaxis()->SetTitle("#bar{#nu}_{e} Pi Estimate / Pi MC");
541  hratioNuePi->GetYaxis()->SetRangeUser(0.94,1.02);
542  hratioNuePi->Draw("hist");
543  line->Draw("same");
544  c1->Write("nue1Dpi");
545  }
546 
547  void FluxDecomp::SavePlots(TDirectory * dir) {
548  if(!fIsDecomposed) Decompose();
549  dir->cd();
550  auto c1 = new TCanvas(UniqueName().c_str(),"c1",800,500);
551  double thePOT=fData.POT();
552 
553  auto hs = new THStack();
554  auto hs2 = new THStack();
555 
556  std::vector<std::pair<Spectrum, int> > sp_col = {
557  {fNueCCFromOther,kGreen+2},
558  {fNueCCFromKaEstim, kMagenta-9},
559  {fNueCCFromPiEstim, kAzure +6}
560  };
561 
562  std::vector<std::pair<Spectrum, int> > sp_col2 = {
563  {fNueCCFromOther,kGreen+2},
564  {fNueCCFromKa, kMagenta -9},
565  {fNueCCFromPi, kAzure+6}
566  };
567 
568  for (auto sp:sp_col){
569  auto hist = sp.first.ToTH1(thePOT, sp.second);
570  hs->Add(hist, "hist");
571  }
572 
573  for (auto sp:sp_col2){
574  auto hist = sp.first.ToTH1(thePOT, sp.second);
575  hist->SetLineStyle(2);
576  hs2->Add(hist, "hist");
577  }
578 
579  //fData.ToTH1(thePOT)->Draw("ep");
580  NueEstimate() .ToTH1(thePOT,kMagenta,1)->Draw("hist");
581  hs->Draw("same");
582  hs2->Draw("same");
583  fAntiNue .ToTH1(thePOT,kMagenta,2)->Draw("hist same");
584  NueEstimate() .ToTH1(thePOT,kMagenta,1)->Draw("hist same");
585  TLegend * leg = new TLegend(0.6, 0.6, 0.89,0.89);
586  TGraph * dummy = new TGraph;
587  dummy->SetLineWidth(3);
588  dummy->SetLineColor(kMagenta); leg->AddEntry(dummy->Clone(),"Total Beam #bar{#nu}_{e}","l");
589  dummy->SetLineColor(kAzure+6); leg->AddEntry(dummy->Clone(),"#bar{#nu}_{e} from Pi","l");
590  dummy->SetLineColor(kMagenta-9); leg->AddEntry(dummy->Clone(),"#bar{#nu}_{e} from Ka","l");
591  dummy->SetLineColor(kGreen+2); leg->AddEntry(dummy->Clone(),"#bar{#nu}_{e} from Other","l");
592  dummy->SetLineColor(kBlack ); dummy->SetLineStyle(2);
593  leg->AddEntry(dummy->Clone(),"Uncorrected MC","l");
594  leg->Draw();
595  c1->Write("Fluxresult");
596  // Nue estimate/mc ratio
597  auto hratioNue = (NueEstimate()/fAntiNue).ToTH1();
598  auto hratioNuePi = (fNueCCFromPiEstim/fNueCCFromPi).ToTH1();
599  auto hratioNueKa = (fNueCCFromKaEstim/fNueCCFromKa).ToTH1();
600  hratioNue->GetYaxis()->SetTitle("#bar{#nu}_{e} Estimate / MC");
601  hratioNue->GetYaxis()->SetRangeUser(0.6,1.2);
602  hratioNue->SetLineColor(kMagenta);
603  hratioNuePi->SetLineColor(kAzure+6);
604  hratioNueKa->SetLineColor(kMagenta-9);
605  hratioNue->Draw("hist");
606  hratioNuePi->Draw("hist same");
607  hratioNueKa->Draw("hist same");
608  leg->GetListOfPrimitives()->RemoveLast();
609  leg->GetListOfPrimitives()->RemoveLast();
610  leg->SetY1(0.15);
611  leg->SetY2(0.45);
612  leg->Draw();
613  c1->Write("Fluxratios");
614  }
615 
616  //----------------------------------------------------------------------
617 
619  return fNue;
620  }
622  return fAntiNue;
623  }
625  return fNumu;
626  }
628  return fAntiNumu;
629  }
631  return fNCTot;
632  }
634  return fNC;
635  }
637  return fNCAnti;
638  }
640  return fData;
641  }
642 }
void NueEstimateFromKa() const
Definition: FluxDecomp.cxx:180
T max(const caf::Proxy< T > &a, T b)
const std::vector< Binning > & GetBinnings() const
Definition: Spectrum.h:256
Spectrum fNumuSelCCFromPi
Definition: FluxDecomp.h:176
const XML_Char * name
Definition: expat.h:151
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
virtual Spectrum MC_AntiNueComponent() const override
Definition: FluxDecomp.cxx:621
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
Spectrum fNC
Definition: FluxDecomp.h:159
Spectrum fNumuSelCCFromOther
Definition: FluxDecomp.h:178
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:148
Spectrum NCAntiComponent() const override
Definition: FluxDecomp.cxx:300
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
Spectrum fNumuSelCCFromKa
Definition: FluxDecomp.h:177
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
Spectrum fTotal
Definition: FluxDecomp.h:166
Spectrum fNumuUncontainBkg
Definition: FluxDecomp.h:182
const Color_t kMC
const Eigen::ArrayXd & GetEigen() const
NB these don&#39;t have POT scaling. For expert high performance ops only!
Definition: Spectrum.h:188
Spectrum fAntiNumu
Definition: FluxDecomp.h:165
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 .
Spectrum fNueCCFromOther
Definition: FluxDecomp.h:171
void abs(TH1 *hist)
virtual Spectrum MC_AntiNumuComponent() const override
Definition: FluxDecomp.cxx:627
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
Spectrum NCTotalComponent() const override
Definition: FluxDecomp.cxx:283
Spectrum fNotNue
Definition: FluxDecomp.h:167
void Scale(double c)
Multiply this spectrum by a constant c.
Definition: Spectrum.cxx:260
virtual Spectrum MC_NumuComponent() const override
Definition: FluxDecomp.cxx:624
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
static std::unique_ptr< FluxDecomp > LoadFrom(TDirectory *dir, const std::string &name)
Definition: FluxDecomp.cxx:383
Spectrum fNumuSelData
Definition: FluxDecomp.h:174
Spectrum AntiNumuComponent() const override
Definition: FluxDecomp.cxx:274
void MakeWeightsNumuFromKaon() const
Definition: FluxDecomp.cxx:192
Spectrum fNumuUncontainCCFromOther
Definition: FluxDecomp.h:185
Spectrum AntiNueComponent() const override
Definition: FluxDecomp.cxx:321
virtual ~FluxDecomp()
Definition: FluxDecomp.cxx:138
Spectrum fNumuUncontainData
Definition: FluxDecomp.h:181
Spectrum Data_Component() const override
Definition: FluxDecomp.cxx:639
double fKaonNormalization
Definition: FluxDecomp.h:192
void NueEstimateFromPi() const
Definition: FluxDecomp.cxx:143
virtual Spectrum MC_NueComponent() const override
Definition: FluxDecomp.cxx:618
Spectrum fNueCCFromPiPtPz
Definition: FluxDecomp.h:172
Spectrum NueComponent() const override
Definition: FluxDecomp.cxx:311
Spectrum fNueCCFromKaEstim
Definition: FluxDecomp.h:188
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:546
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
virtual Spectrum MC_NCAntiComponent() const override
Definition: FluxDecomp.cxx:636
void SavePlotsPi(TDirectory *dir)
Definition: FluxDecomp.cxx:512
Spectrum fNueCCFromPi
Definition: FluxDecomp.h:169
void SavePlots(TDirectory *dir)
Definition: FluxDecomp.cxx:547
virtual Spectrum MC_NCTotalComponent() const override
Definition: FluxDecomp.cxx:630
#define pot
Spectrum fNue
Definition: FluxDecomp.h:162
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:517
Spectrum fNumuUncontainCCFromPi
Definition: FluxDecomp.h:183
void SavePlotsKa(TDirectory *dir)
Definition: FluxDecomp.cxx:429
std::vector< float > Spectrum
Definition: Constants.h:573
double POT() const
Definition: Spectrum.h:219
OStream cout
Definition: OStream.cxx:6
Spectrum fNumuSelCCFromPiPtPz
Definition: FluxDecomp.h:179
Base class for the various types of spectrum loader.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const Cut cut
Definition: exporter_fd.C:30
void Decompose() const
Definition: FluxDecomp.cxx:129
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
TH3 * ToTH3(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 3D to obtain TH3.
Definition: Spectrum.cxx:218
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
Spectrum fNCAnti
Definition: FluxDecomp.h:160
const Cut kIsAllFromPions([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return( (sr->mc.nu[0].beam.ndecay==11|| sr->mc.nu[0].beam.ndecay==12|| sr->mc.nu[0].beam.ndecay==13|| sr->mc.nu[0].beam.ndecay==14)&& ((abs(sr->mc.nu[0].beam.tptype)==211)|| (abs(sr->mc.nu[0].beam.tptype)==3122)|| (abs(sr->mc.nu[0].beam.tptype)==3222)|| (abs(sr->mc.nu[0].beam.tptype)==3112)|| (abs(sr->mc.nu[0].beam.tptype)==3212)) );})
Definition: BeamCuts.h:35
Int_t nbinsx
Definition: plot.C:23
TDirectory * dir
Definition: macro.C:5
REGISTER_LOADFROM("BENDecomp", IDecomp, BENDecomp)
Spectrum fNumuSelBkg
Definition: FluxDecomp.h:175
virtual Spectrum MC_NCComponent() const override
Definition: FluxDecomp.cxx:633
Spectrum fNueCCFromKa
Definition: FluxDecomp.h:170
double epsilon
std::vector< Loaders * > loaders
Definition: syst_header.h:386
float fNC[pidBins]
Spectrum NumuComponent() const override
Definition: FluxDecomp.cxx:266
TH2 * fNumuPiWeights
Definition: FluxDecomp.h:191
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
c1
Definition: demo5.py:24
Spectrum fNumuUncontainCCFromKa
Definition: FluxDecomp.h:184
Spectrum fNCTot
Definition: FluxDecomp.h:158
This module creates Common Analysis Files.
Definition: FileReducer.h:10
Spectrum fNueCCFromPiEstim
Definition: FluxDecomp.h:187
T min(const caf::Proxy< T > &a, T b)
Spectrum fAntiNue
Definition: FluxDecomp.h:163
const std::vector< std::string > & GetLabels() const
Definition: Spectrum.h:255
Double_t sum
Definition: plot.C:31
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
void MakeWeightsNumuFromPion() const
Definition: FluxDecomp.cxx:236
virtual Spectrum NueEstimate() const
Definition: FluxDecomp.cxx:111
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
Spectrum NCComponent() const override
Definition: FluxDecomp.cxx:291
Spectrum fData
Definition: FluxDecomp.h:157
void SaveTo(TDirectory *dir, const std::string &name) const override
Definition: FluxDecomp.cxx:330
Spectrum fNumu
Definition: FluxDecomp.h:164