nueBkgEMComposition.C
Go to the documentation of this file.
1 #include "OscLib/OscCalcDumb.h"
2 
3 #include "CAFAna/Core/Spectrum.h"
7 
8 #include "CAFAna/Cuts/Cuts.h"
12 
13 #include "CAFAna/Systs/XSecSysts.h"
14 
15 #include "CAFAna/Vars/Vars.h"
20 
23 
25 
26 #include "Utilities/func/MathUtil.h"
27 
28 #include "TCanvas.h"
29 #include "TH2.h"
30 #include "TH1.h"
31 #include "TFile.h"
32 #include "TLegend.h"
33 #include "TLorentzVector.h"
34 #include "TArrow.h"
35 #include "TMath.h"
36 #include "TLatex.h"
37 #include "TPolyLine.h"
38 
39 #include <iostream>
40 #include <cmath>
41 
42 using namespace ana;
43 
44 // Need a custom var to check whether we have a prim pi0 in event
45 // If we do, only count it if it has > 500 MeV of momentum
46 const Var kNPi0([](const caf::SRProxy* sr){
47  int npi = 0;
48  if (sr->mc.nnu == 0) return npi;
49  for (auto &p : sr->mc.nu[0].prim){
50  if (p.pdg == 111 && p.p.T() > sqrt(0.5*0.5+0.135*0.135)) npi++;
51  }
52  return npi;
53  });
54 
55 
57 {
58  TLatex* prelim = new TLatex(.9, .95, "NOvA Simulation");
59  prelim->SetTextColor(kGray+1);
60  prelim->SetNDC();
61  prelim->SetTextSize(2/30.);
62  prelim->SetTextAlign(32);
63  prelim->Draw();
64 }
65 
66 std::vector<double> MakeTable(bool isRHC, bool is0pi)
67 {
68  std::string fname = "Ana2018_bkgd_pi0cut.root";
69  std::string appen = isRHC ? "RHC" : "FHC";
70  if (is0pi) appen += "_0pi";
71  else appen += "_pi";
72 
73  // Grab predictions from our file
74  IPrediction *pred = LoadFromFile<IPrediction>(fname,"pred"+appen).release();
75 
76  // Grab our exposures
77  double pot = isRHC ? kAna2018RHCPOT : kAna2018FHCPOT;
78  double lt = isRHC ? kAna2018RHCLivetime : kAna2018FHCLivetime;
79 
80  // Set up oscillations, set to 2018 best fit
81  double dCP = 0.166;
82  double ssth23 = 0.584;
83  double ssth23UO = ssth23;
84  double ssth23LO = 0.514 - (ssth23-0.514);
85  double dmsq32 = 0.002505;
86 
90 
91  kFitDeltaInPiUnits.SetValue(calc,dCP);
92  kFitDmSq32.SetValue(calc,dmsq32);
93  kFitSinSqTheta23.SetValue(calc,ssth23);
94 
95  double wsBF = pred->PredictComponent(calc, ana::Flavors::kNuMuToNuE,
98  if (!isRHC)
99  wsBF = pred->PredictComponent(calc, ana::Flavors::kNuMuToNuE,
102 
103  if (isRHC) {
104  // High = NH, 3pi/2 UO
105  // Low = IH, pi/2, LO
106  kFitDeltaInPiUnits.SetValue(calcL,0.5);
107  kFitDeltaInPiUnits.SetValue(calcH,1.5);
108  kFitSinSqTheta23.SetValue(calcL,ssth23LO);
109  kFitSinSqTheta23.SetValue(calcH,ssth23UO);
110  kFitDmSq32.SetValue(calcL,-dmsq32);
111  kFitDmSq32.SetValue(calcH,dmsq32);
112  }
113  else {
114  // High = IH, pi/2, UO
115  // Low = NH, 3pi/2, LO
116  kFitDeltaInPiUnits.SetValue(calcL,1.5);
117  kFitDeltaInPiUnits.SetValue(calcH,0.5);
118  kFitSinSqTheta23.SetValue(calcL,ssth23LO);
119  kFitSinSqTheta23.SetValue(calcH,ssth23UO);
120  kFitDmSq32.SetValue(calcL,dmsq32);
121  kFitDmSq32.SetValue(calcH,-dmsq32);
122  }
123 
124 
125  // Make predictions for each background component
126  TH1 *gne = pred->PredictComponent(calc, ana::Flavors::kNuEToNuE,
128  ana::Sign::kNu).ToTH1(pot);
129  TH1 *gnebar = pred->PredictComponent(calc, ana::Flavors::kNuEToNuE,
132  // Merge beam nue components
133  gne->Add(gnebar);
134 
135 
136  TH1 *gwsL = pred->PredictComponent(calcL, ana::Flavors::kNuMuToNuE,
139  ).ToTH1(pot);
140  TH1 *gwsH = pred->PredictComponent(calcH, ana::Flavors::kNuMuToNuE,
142  (isRHC ? ana::Sign::kNu : ana::Sign::kAntiNu)
143  ).ToTH1(pot);
144 
145  TH1 *gnm = pred->PredictComponent(calc, ana::Flavors::kNuMuToNuMu,
147  ana::Sign::kBoth).ToTH1(pot);
148  TH1 *gnc = pred->PredictComponent(calc, ana::Flavors::kAll,
150  ana::Sign::kBoth).ToTH1(pot);
151  TH1 *gtau= pred->PredictComponent(calc, ana::Flavors::kNuMuToNuTau,
153  ana::Sign::kBoth).ToTH1(pot);
154  gtau->Add(pred->PredictComponent(calc, ana::Flavors::kNuEToNuTau,
156  ana::Sign::kBoth).ToTH1(pot));
157 
158  std::cout << "WS range: " << gwsL->Integral() << " to " << gwsH->Integral() << std::endl;
159  std::vector<double> vals = {gnm->Integral(),gnc->Integral(),
160  gne->Integral(),gwsL->Integral(),
161  gwsH->Integral(),
162  gtau->Integral(),
163  wsBF};
164  return vals;
165 }
166 
168 {
169  std::string appen = isRHC ? "RHC" : "FHC";
170 
171  std::vector<double> vals_0pi = MakeTable(isRHC, true);
172  std::vector<double> vals_pi = MakeTable(isRHC, false);
173 
174  TCanvas *can = new TCanvas(UniqueName().c_str(),"",1200,900);
175 
176  double topmargin = 0.1;
177  double botmargin = 0.12;
178  gPad->SetBottomMargin(botmargin);
179  gPad->SetBottomMargin(topmargin);
180 
181  double xmax = isRHC ? 5 : 9;
182  TH2 *axes = new TH2F(UniqueName().c_str(),"",1,0,xmax,1,0,4);
183  axes->GetXaxis()->SetTitle("Number of Bkg Events");
184  axes->GetXaxis()->CenterTitle();
185  axes->GetXaxis()->SetTitleSize(0.06);
186  axes->GetXaxis()->SetTitleOffset(0.7);
187  axes->GetYaxis()->SetBinLabel(1,"");
188  axes->Draw();
189 
190  TLatex *l1 = new TLatex(0.07,0.33,"#nu_{e} CC");
191  l1->SetNDC();
192  l1->SetTextSize(0.06);
193  l1->SetTextAngle(90);
194  l1->Draw();
195  TLatex *l2 = new TLatex(0.07,0.56, "NC");
196  l2->SetNDC();
197  l2->SetTextSize(0.06);
198  l2->SetTextAngle(90);
199  l2->Draw();
200  TLatex *l3 = new TLatex(0.07,0.725,"#nu_{#mu} CC");
201  l3->SetNDC();
202  l3->SetTextSize(0.06);
203  l3->SetTextAngle(90);
204  l3->Draw();
205  TLatex *l4 = new TLatex(0.07,0.105,"Cosmic");
206  l4->SetNDC();
207  l4->SetTextSize(0.06);
208  l4->SetTextAngle(90);
209  l4->Draw();
210 
211 
212  std::string potStr =
213  isRHC ? "6.91#times10^{20} POT-equiv" : "8.85#times10^{20} POT-equiv";
214  TLatex *lpot = new TLatex(0.45,0.80,potStr.c_str());
215  lpot->SetNDC();
216  lpot->SetTextSize(0.06);
217  lpot->Draw();
218 
219 
220 
221  double buffer = 0.15;
222 
223  TBox *nm0pi = new TBox(0,3+buffer,vals_0pi[0],4-buffer);
224  TBox *nmpi = new TBox(vals_0pi[0],3+buffer,vals_0pi[0]+vals_pi[0],4-buffer);
225  nm0pi->SetFillColor(kOrange-3);
226  nmpi->SetFillColor(kAzure+1);
227  nm0pi->Draw("same");
228  nmpi->Draw("same");
229 
230  TBox *nc0pi = new TBox(0,2+buffer,vals_0pi[1],3-buffer);
231  TBox *ncpi = new TBox(vals_0pi[1],2+buffer,vals_0pi[1]+vals_pi[1],3-buffer);
232  nc0pi->SetFillColor(kOrange-3);
233  ncpi->SetFillColor(kAzure+1);
234  nc0pi->Draw("same");
235  ncpi->Draw("same");
236 
237  double lowE;
238  double highE;
239  lowE = 0;
240  highE = lowE + vals_0pi[2] + vals_pi[2];
241  TBox *nelec = new TBox(lowE,1+buffer,highE,2-buffer);
242 
243  lowE = highE;
244  highE = lowE + vals_0pi[3] + vals_pi[3];
245  TBox *nws = new TBox(lowE,1+buffer,highE,2-buffer);
246 
247  lowE = highE;
248  highE = lowE + vals_0pi[4] + vals_pi[4];
249  TBox *nwsH = new TBox(lowE-0.02,1+buffer+0.013,highE,2-buffer-0.013);
250 
251 
252  nelec->SetFillColor(kMagenta);
253  if (isRHC){
254  nws->SetFillColor(kViolet-5); // nue app color
255  nwsH->SetLineColor(kViolet-5); // nue app color
256  }
257  else{
258  nws->SetFillColor(kMagenta+3); // nuebar app color
259  nwsH->SetLineColor(kMagenta+3); // nuebar app color
260  }
261  nwsH->SetLineStyle(7);
262  nwsH->SetLineWidth(3);
263  nwsH->SetFillStyle(0);
264 
265  nelec->Draw("same");
266  nwsH->Draw("F");
267  nws->Draw("same");
268 
269  // Making the arrow
270  double arrowBuffer = isRHC ? 0.2*5/9 : 0.2;
271  Double_t triLX[4] = {lowE,lowE+arrowBuffer,lowE+arrowBuffer,lowE};
272  Double_t triLY[4] = {1.5,1.6,1.4,1.5};
273  TPolyLine *triL = new TPolyLine(4, triLX, triLY);
274 
275  Double_t triRX[4] = {highE,highE-arrowBuffer,highE-arrowBuffer,highE};
276  Double_t triRY[4] = {1.5,1.6,1.4,1.5};
277  TPolyLine *triR = new TPolyLine(4, triRX, triRY);
278 
279  TLine *lArrow = new TLine(lowE+0.1,1.5,highE-0.1,1.5);
280  lArrow->SetLineWidth(3);
281 
282  if (isRHC){
283  triR->SetFillColor(kViolet-5); // nue app color
284  triL->SetFillColor(kViolet-5); // nue app color
285  lArrow->SetLineColor(kViolet-5);// nue app color
286  }
287  else{
288  triR->SetFillColor(kMagenta+3); // nue app color
289  triL->SetFillColor(kMagenta+3); // nue app color
290  lArrow->SetLineColor(kMagenta+3);// nue app color
291  }
292  triL->SetLineColor(0);
293  triR->SetLineColor(0);
294  triL->SetFillStyle(1001);
295  triR->SetFillStyle(1001);
296 
297  triR->Draw("f");
298  triL->Draw("f");
299  lArrow->Draw("same");
300 
301  std::string beam = isRHC ? "rhc" : "fhc";
302  std::pair <TH1D*, double > hcosmic = GetNueCosmics2018(beam);
303  double cosmic = hcosmic.first->Integral();
304 
305  TBox *cos = new TBox(0,0+buffer,cosmic,1-buffer);
306  cos->SetFillColor(kGray+2);
307  cos->Draw("same");
308 
309 
310  TLegend *leg = new TLegend(0.40,0.50,0.55,0.78);
311  leg->AddEntry(nelec,"Beam #nu_{e} and #bar{#nu}_{e}","f");
312  if (isRHC)
313  leg->AddEntry(nws,"Wrong-sign (Appearing #nu_{e})","f");
314  else
315  leg->AddEntry(nws,"Wrong-sign (Appearing #bar{#nu}_{e})","f");
316  leg->AddEntry(nmpi, "Beam bkg with a #pi^{0} (KE > 0.5 GeV)","f");
317  leg->AddEntry(nm0pi,"Beam bkg without EM activity","f");
318  leg->AddEntry(cos, "Cosmic bkg","f");
319  leg->SetTextSize(0.04);
320  leg->SetFillStyle(0);
321  leg->Draw();
322 
323  gPad->RedrawAxis();
324 
325  PrintSimulation();
326 
327  std::string modeText = isRHC ? "Antineutrino Beam" : "Neutrino Beam";
328  TLatex *text = new TLatex(0.1,0.92,modeText.c_str());
329  text->SetNDC();
330  text->SetTextSize(0.04);
331  text->Draw();
332 
333  can->SaveAs(("Fig_FDBkgd_Pi0Comp_"+appen+".pdf").c_str());
334 
335 }
336 
338 {
339  SpectrumLoader lSwapRHC("prod_sumdecaf_R17-11-14-prod4reco.e_fd_genie_fluxswap_rhc_nova_v08_full_v1_nue2018");
340  SpectrumLoader lNonSRHC("prod_sumdecaf_R17-11-14-prod4reco.e_fd_genie_nonswap_rhc_nova_v08_full_v1_nue2018");
341  SpectrumLoader lTauSRHC("prod_sumdecaf_R17-11-14-prod4reco.e_fd_genie_tau_rhc_nova_v08_full_v1_nue2018");
342 
346 
350 
351  PredictionNoExtrap predRHC_0pi(lNonSRHC,lSwapRHC,lTauSRHC,"",kNue2018Binning,
352  kNue2018AnaBin,sel&&kNPi0==0,kNoShift,wei);
353  PredictionNoExtrap predRHC_pi (lNonSRHC,lSwapRHC,lTauSRHC,"",kNue2018Binning,
354  kNue2018AnaBin,sel&&kNPi0>0,kNoShift,wei);
355 
356 
357  SpectrumLoader lSwapFHC("prod_sumdecaf_R17-11-14-prod4reco.d_fd_genie_fluxswap_fhc_nova_v08_full_v1_nue2018");
358  SpectrumLoader lNonSFHC("prod_sumdecaf_R17-11-14-prod4reco.d_fd_genie_nonswap_fhc_nova_v08_full_v1_nue2018");
359  SpectrumLoader lTauSFHC("prod_sumdecaf_R17-11-14-prod4reco.d_fd_genie_tau_fhc_nova_v08_full_v1_nue2018");
360 
364 
365  PredictionNoExtrap predFHC_0pi(lNonSFHC,lSwapFHC,lTauSFHC,"",kNue2018Binning,
366  kNue2018AnaBin,sel&&kNPi0==0,kNoShift,wei);
367  PredictionNoExtrap predFHC_pi (lNonSFHC,lSwapFHC,lTauSFHC,"",kNue2018Binning,
368  kNue2018AnaBin,sel&&kNPi0>0,kNoShift,wei);
369 
370  lNonSRHC.Go();
371  lTauSRHC.Go();
372  lSwapRHC.Go();
373 
374  lNonSFHC.Go();
375  lTauSFHC.Go();
376  lSwapFHC.Go();
377 
378 
379  TFile* outFile = new TFile("Ana2018_bkgd_pi0cut.root","RECREATE");
380  outFile->cd();
381 
382  predFHC_0pi.SaveTo(outFile, "predFHC_0pi");
383  predFHC_pi. SaveTo(outFile, "predFHC_pi");
384 
385  predRHC_0pi.SaveTo(outFile, "predRHC_0pi");
386  predRHC_pi. SaveTo(outFile, "predRHC_pi");
387 
388  outFile->Close();
389 }
390 
391 
392 void nueBkgEMComposition(bool fillSpectra)
393 {
394  if (fillSpectra){
395  FillSpectra();
396  }
397  else{
398  MakeBkgEMPlot(true);
399  MakeBkgEMPlot(false);
400  }
401 }
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::map< std::string, double > xmax
virtual Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
double lt
double ssth23
Antineutrinos-only.
Definition: IPrediction.h:50
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
(&#39; appearance&#39;)
Definition: IPrediction.h:18
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
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
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:18
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Definition: FitVars.cxx:177
const Var kNPi0([](const caf::SRProxy *sr){int npi=0;if(sr->mc.nnu==0) return npi;for(auto &p:sr->mc.nu[0].prim){if(p.pdg==111 &&p.p.T() > sqrt(0.5 *0.5+0.135 *0.135)) npi++;}return npi;})
(&#39;beam &#39;)
Definition: IPrediction.h:15
const char * p
Definition: xmltok.h:285
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2043
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:570
caf::Proxy< short int > nnu
Definition: SRProxy.h:569
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Forward to wrapped Var&#39;s SetValue()
void SetSpillCut(const SpillCut &cut)
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:333
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
::xsd::cxx::tree::buffer< char > buffer
Definition: Database.h:179
const double kAna2018RHCPOT
Definition: Exposures.h:208
const Var kNue2018AnaBin([](const caf::SRProxy *sr){int selBin=kNue2018SelectionBin(sr);float nuE=kNueEnergy2018(sr);int nuEBin=nuE/0.5;assert(nuEBin<=8 &&"An event with nuE > 4.5 should never happen");int anaBin=9 *selBin+nuEBin;return anaBin;})
Use this Analysis Binning for Ana2018, official Binning.
Definition: NueCuts2018.h:171
Charged-current interactions.
Definition: IPrediction.h:39
Sum up livetimes from individual cosmic triggers.
TFile * outFile
Definition: PlotXSec.C:135
void MakeBkgEMPlot(bool isRHC)
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
double dCP
#define pot
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
const Binning kNue2018Binning
Definition: NueCuts2018.cxx:92
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
TLatex * prelim
Definition: Xsec_final.C:133
Neutrinos-only.
Definition: IPrediction.h:49
const SystShifts kNoShift
Definition: SystShifts.h:115
OStream cout
Definition: OStream.cxx:6
(&#39; survival&#39;)
Definition: IPrediction.h:19
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2055
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::pair< TH1D *, double > GetNueCosmics2018(std::string beam, bool GetFromUPS=false, bool NERSC=false)
double dmsq32
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
void PrintSimulation()
const Var kXSecCVWgt2018
Definition: XsecTunes.h:49
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
const Cut kNue2018FDAllSamples
Definition: NueCuts2018.h:77
Neutral-current interactions.
Definition: IPrediction.h:40
T cos(T number)
Definition: d0nt_math.hpp:78
std::vector< double > MakeTable(bool isRHC, bool is0pi)
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:46
void nueBkgEMComposition(bool fillSpectra)
Prediction that just uses FD MC, with no extrapolation.
const double kAna2018FHCPOT
Definition: Exposures.h:207
void FillSpectra()
All neutrinos, any flavor.
Definition: IPrediction.h:26
#define for
Definition: msvc_pragmas.h:3
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
const HistAxis kNue2018Axis("NuE Energy / Analysis Bin", kNue2018Binning, kNue2018AnaBin)
Use this Axis for Ana2018, official Axis.
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:30
const double kAna2018FHCLivetime
Definition: Exposures.h:213
const double kAna2018RHCLivetime
Definition: Exposures.h:214