CompareMissingLeptons.C
Go to the documentation of this file.
1 /// CompareMissingLeptons.C by J. Hewes <jhewes15@fnal.gov>
2 /// Compare NC spectra with and without a missing lepton cut in MC
3 
5 
8 
9 #include "TFile.h"
10 #include "TCanvas.h"
11 
12 using namespace ana;
13 
14 using std::cout;
15 using std::endl;
16 using std::string;
17 using std::vector;
18 
19 void DrawSpectrum(TH1D* h1, TH1D* h2, string name, vector<double> errors={}) {
20 
21  cout << "With leptons integral (" << name << "): " << h1->Integral() << endl;
22  cout << "No leptons integral (" << name << "): " << h2->Integral() << endl;
23 
24  TCanvas c("c", "c", 1600, 900);
25  h1->SetLineWidth(2);
26  h1->SetLineColor(kRed-3);
27  h2->SetLineWidth(2);
28  h2->SetLineColor(kBlue-3);
29 
30  h1->Draw("hist");
31  h2->Draw("hist same");
32 
33  string plotName = "spectrum_" + name + ".pdf";
34  c.SaveAs(plotName.c_str());
35 
36  h2->Divide(h1);
37  h2->SetLineColor(kBlack);
38 
39  double max = fabs(1-h2->GetMinimum());
40  max = 0.1 * ceil(max * 10);
41  h2->SetMinimum(1-max);
42  h2->SetMaximum(1);
43 
44  if (errors.size() > 0) {
45  for (size_t i = 0; i < errors.size(); ++i) {
46  h2->SetBinError(i+1, errors[i] * h2->GetBinContent(i+1));
47  }
48  h2->Draw("hist");
49  } else {
50  h2->Draw("hist");
51  }
52 
53  plotName = "ratio_" + name + ".pdf";
54  c.SaveAs(plotName.c_str());
55 
56 }
57 
59 
60  IPrediction* predWithLeptons = LoadFromFile<FDPredictionSterile>("/pnfs/nova/persistent/analysis/nux/binning_studies/pred_ncsel_fhc_fardet.root", "pred_ncsel_fhc_fardet").release();
61  IPrediction* predNoLeptons = LoadFromFile<FDPredictionSterile>("/pnfs/nova/persistent/analysis/nux/binning_studies/pred_ncsel_fhc_fardet_missingleptoncut.root", "pred_ncsel_fhc_fardet").release();
62 
63  auto calc = DefaultSterileCalc(4);
65  // calc->SetDm(4, 10);
66  // calc->SetAngle(2, 4, 0.785398);
67  // calc->SetAngle(3, 4, 1.5708);
68 
69  TH1D* hWithLeptons = predWithLeptons->Predict(calc).ToTH1(9.5e20);
70  TH1D* hNoLeptons = predNoLeptons->Predict(calc).ToTH1(9.5e20);
71  DrawSpectrum(hWithLeptons, hNoLeptons, "full");
72 
73  hWithLeptons = predWithLeptons->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(9.5e20);
74  hNoLeptons = predNoLeptons ->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(9.5e20);
77  TH1D* hWL = sWL.ToTH1(sWL.POT());
78  TH1D* hNL = sNL.ToTH1(sNL.POT());
79  vector<double> errors;
80  for (size_t i = 1; i <= hWL->GetNbinsX(); ++i) {
81  errors.push_back(sqrt((1/hWL->GetBinContent(i)) + (1/hNL->GetBinContent(i))));
82  }
83  DrawSpectrum(hWithLeptons, hNoLeptons, "signal", errors);
84 
85  hWithLeptons = predWithLeptons->PredictComponent(calc, Flavors::kAll, Current::kCC, Sign::kBoth).ToTH1(9.5e20);
86  hNoLeptons = predNoLeptons ->PredictComponent(calc, Flavors::kAll, Current::kCC, Sign::kBoth).ToTH1(9.5e20);
87  DrawSpectrum(hWithLeptons, hNoLeptons, "background");
88 
89  hWithLeptons = predWithLeptons->PredictComponent(calc, Flavors::kNuEToNuE, Current::kCC, Sign::kBoth).ToTH1(9.5e20);
90  hNoLeptons = predNoLeptons ->PredictComponent(calc, Flavors::kNuEToNuE, Current::kCC, Sign::kBoth).ToTH1(9.5e20);
91  DrawSpectrum(hWithLeptons, hNoLeptons, "nue_to_nue");
92 
93  hWithLeptons = predWithLeptons->PredictComponent(calc, Flavors::kNuMuToNuMu, Current::kCC, Sign::kBoth).ToTH1(9.5e20);
94  hNoLeptons = predNoLeptons ->PredictComponent(calc, Flavors::kNuMuToNuMu, Current::kCC, Sign::kBoth).ToTH1(9.5e20);
95  DrawSpectrum(hWithLeptons, hNoLeptons, "numu_to_numu");
96 
97  hWithLeptons = predWithLeptons->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kBoth).ToTH1(9.5e20);
98  hNoLeptons = predNoLeptons ->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kBoth).ToTH1(9.5e20);
99  DrawSpectrum(hWithLeptons, hNoLeptons, "numu_to_nue");
100 
101  hWithLeptons = predWithLeptons->PredictComponent(calc, Flavors::kNuEToNuMu, Current::kCC, Sign::kBoth).ToTH1(9.5e20);
102  hNoLeptons = predNoLeptons ->PredictComponent(calc, Flavors::kNuEToNuMu, Current::kCC, Sign::kBoth).ToTH1(9.5e20);
103  DrawSpectrum(hWithLeptons, hNoLeptons, "nue_to_numu");
104 
105  hWithLeptons = predWithLeptons->PredictComponent(calc, Flavors::kNuMuToNuTau, Current::kCC, Sign::kBoth).ToTH1(9.5e20);
106  hNoLeptons = predNoLeptons ->PredictComponent(calc, Flavors::kNuMuToNuTau, Current::kCC, Sign::kBoth).ToTH1(9.5e20);
107  DrawSpectrum(hWithLeptons, hNoLeptons, "numu_to_nutau");
108 
109  hWithLeptons = predWithLeptons->PredictComponent(calc, Flavors::kNuEToNuTau, Current::kCC, Sign::kBoth).ToTH1(9.5e20);
110  hNoLeptons = predNoLeptons ->PredictComponent(calc, Flavors::kNuEToNuTau, Current::kCC, Sign::kBoth).ToTH1(9.5e20);
111  DrawSpectrum(hWithLeptons, hNoLeptons, "nue_to_nutau");
112 
113 } // macro CompareMissingLeptons
const XML_Char * name
Definition: expat.h:151
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
virtual Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp: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:148
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
(&#39;beam &#39;)
Definition: IPrediction.h:15
T sqrt(T number)
Definition: d0nt_math.hpp:156
void DrawSpectrum(TH1D *h1, TH1D *h2, string name, vector< double > errors={})
osc::OscCalcDumb calc
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
Charged-current interactions.
Definition: IPrediction.h:39
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
void SetNus20Params(osc::OscCalcSterile *calc)
Definition: Utilities.h:682
double POT() const
Definition: Spectrum.h:227
OStream cout
Definition: OStream.cxx:6
(&#39; survival&#39;)
Definition: IPrediction.h:19
void CompareMissingLeptons()
TH1F * h2
Definition: plot.C:45
TH1F * h1
Neutral-current interactions.
Definition: IPrediction.h:40
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
All neutrinos, any flavor.
Definition: IPrediction.h:26
(&#39; appearance&#39;)
Definition: IPrediction.h:16
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
enum BeamMode kBlue
fvar< T > ceil(const fvar< T > &x)
Definition: ceil.hpp:11
enum BeamMode string