numu_sig_nonmax.C
Go to the documentation of this file.
2 #include "OscLib/OscCalcPMNS.h"
4 #include "OscLib/OscCalc.h"
5 
6 #include "CAFAna/Cuts/Cuts.h"
9 #include "CAFAna/Vars/FitVars.h"
11 #include "CAFAna/Analysis/Plots.h"
15 #include "CAFAna/Core/Utilities.h"
16 #include "CAFAna/Vars/Vars.h"
18 #include "CAFAna/Core/Binning.h"
20 
21 #include "TCanvas.h"
22 #include "TFile.h"
23 #include "TGraph.h"
24 #include "TH1.h"
25 #include "TH2.h"
26 #include "TMath.h"
27 #include "TGaxis.h"
28 #include "TMultiGraph.h"
29 #include "TLegend.h"
30 #include "TLatex.h"
31 #include "TStyle.h"
32 #include "THStack.h"
33 #include "TPaveText.h"
34 
35 #include <cmath>
36 #include <iostream>
37 #include <sstream>
38 #include <string>
39 #include <fstream>
40 #include <iomanip>
41 
42 #include "Utilities/func/MathUtil.h"
43 
44 // No Cosmics, no Systematics, S15-05-22 hadded cafs
45 // Use with Nova style ROOT logon script
46 
47 using namespace ana;
48 
49 void Simulation() // for some reason doesn't always remember definition in Style script
50 {
51  TLatex* prelim = new TLatex(.9, .95, "NO#nuA Simulation");
52  prelim->SetTextColor(kGray+1);
53  prelim->SetNDC();
54  prelim->SetTextSize(2/30.);
55  prelim->SetTextAlign(32);
56  prelim->Draw();
57 }
58 
60 {
61 
62  TGaxis::SetMaxDigits(4);
63 
64  gStyle->SetTitleOffset(.7, "y");
65  gStyle->SetTitleOffset(.77, "x");
66  gStyle->SetOptTitle(0);
67 
68  // Caltech
69  // const std::string fname = "$NOVA_DATA/mc/forFA/fdgenie.nonswap.decaf.root";
70  // const std::string fnameSwap = "$NOVA_DATA/mc/forFA/fdgenie.swap.decaf.root";
71 
72  // FNAL
73  const std::string fname = "/nova/ana/users/bays/files/decafs/fdgenie.nonswap.decaf.root";
74  const std::string fnameSwap = "/nova/ana/users/bays/files/decafs/fdgenie.swap.decaf.root";
75 
77 
78  calc.SetL(810);
79  calc.SetRho(2.75);
80  calc.SetDmsq21(7.5e-5);
81  calc.SetDmsq32(2.35e-3);
82  calc.SetTh12(.601);
83  calc.SetTh13(.1567);
84  calc.SetdCP(0);
85  calc.SetTh23(TMath::Pi()/4); // sin^2(2theta) = 1
86 
87  SpectrumLoader loader(fname);
88  SpectrumLoader loaderSwap(fnameSwap);
89 
90  // Load up the necessary histograms, seperated by flavour etc
91  PredictionNoExtrap predCC(loader, loaderSwap,
92  "Non-QE Neutrino Energy (GeV)", kNumuEnergyBinning,
93  kCCE, kNumuFD);
94 
95  loader.Go();
96  loaderSwap.Go();
97 
98  double temp = 0.0, temp2 = 0.0;;
99 
100  double chisqmap[30] = {0.};
101  double chisqmap2[30] = {0.};
102  double chisqmap3[30] = {0.};
103 
104  double pot = 18e20;
105 
106  for(unsigned int i = 0; i < 30; i++) {
107  temp2 = 0.355+(double)(i)/100.0;
108  temp = asin(sqrt(temp2));
109  calc.SetTh23(temp); // 0.35 to 0.65
110  calc.SetDmsq32(2.35e-3);
111 
112  std::cout << "Filling fake data for : " << temp << " " << temp2 << std::endl;
113 
114  Spectrum obsCC = predCC.Predict(&calc);
115  Spectrum fakeCC = obsCC.FakeData(pot);
116 
117  SingleSampleExperiment exptCC(&predCC, fakeCC);
118 
119  calc.SetTh23(TMath::Pi()/4);
120  MinuitFitter fit(&exptCC,{&kFitDmSq32});
121  chisqmap[i]=fit.Fit(&calc)->EvalMetricVal();
122  }
123 
124  pot = 6e20;
125 
126  for(unsigned int i = 0; i < 30; i++) {
127  temp2 = 0.355+(double)(i)/100.0;
128  temp = asin(sqrt(temp2));
129  calc.SetTh23(temp); // 0.35 to 0.65
130  calc.SetDmsq32(2.35e-3);
131 
132  std::cout << "Filling fake data for : " << temp << " " << temp2 << std::endl;
133 
134  Spectrum obsCC = predCC.Predict(&calc);
135  Spectrum fakeCC = obsCC.FakeData(pot);
136 
137  SingleSampleExperiment exptCC(&predCC, fakeCC);
138 
139  calc.SetTh23(TMath::Pi()/4);
140  MinuitFitter fit(&exptCC,{&kFitDmSq32});
141  chisqmap2[i]=fit.Fit(&calc)->EvalMetricVal();
142  }
143 
144  pot = 2e20;
145 
146  for(unsigned int i = 0; i < 30; i++) {
147  temp2 = 0.355+(double)(i)/100.0;
148  temp = asin(sqrt(temp2));
149  calc.SetTh23(temp); // 0.35 to 0.65
150  calc.SetDmsq32(2.35e-3);
151 
152  std::cout << "Filling fake data for : " << temp << " " << temp2 << std::endl;
153 
154  Spectrum obsCC = predCC.Predict(&calc);
155  Spectrum fakeCC = obsCC.FakeData(pot);
156 
157  SingleSampleExperiment exptCC(&predCC, fakeCC);
158 
159  calc.SetTh23(TMath::Pi()/4);
160  MinuitFitter fit(&exptCC,{&kFitDmSq32});
161  chisqmap3[i]=fit.Fit(&calc)->EvalMetricVal();
162  }
163 
164  TH1F *chisqdiff = new TH1F("chisqdiff","chisqdiff",30,0.35,0.65);
165  TH1F *signif = new TH1F("signif","signif",30,0.35,0.65);
166  TH1F *signif2 = new TH1F("signif2","signif2",30,0.35,0.65);
167  TH1F *signif3 = new TH1F("signif3","signif3",30,0.35,0.65);
168 
169  for (unsigned int i = 1; i<=30; i++) {
170  chisqdiff->SetBinContent(i,chisqmap[i-1]);
171  signif->SetBinContent(i,sqrt(chisqmap[i-1]));
172  signif2->SetBinContent(i,sqrt(chisqmap2[i-1]));
173  signif3->SetBinContent(i,sqrt(chisqmap3[i-1]));
174  }
175 
176  chisqdiff->SetLineWidth(2);
177  chisqdiff->SetLineColor(2);
178 
179  signif->SetLineWidth(2);
180  signif->SetLineColor(4);
181  signif->GetYaxis()->SetTitle("Significance (#sigma)");
182  signif->GetXaxis()->SetTitle("True sin^{2}(#theta_{23})");
183  signif->GetYaxis()->CenterTitle();
184  signif->GetZaxis()->CenterTitle();
185 
186  signif2->SetLineWidth(2);
187  signif2->SetLineColor(2);
188 
189  signif3->SetLineWidth(2);
190  signif3->SetLineColor(1);
191 
192  TLegend *leg = new TLegend(0.25,0.46,0.98,0.88);
193  leg->SetFillStyle(0);
194  leg->SetBorderSize(0);
195  leg->SetHeader("14 kton, signif. to rule out max. mixing");
196 
197  signif->Draw("l");
198  signif2->Draw("lsame");
199  signif3->Draw("lsame");
200 
201  leg->AddEntry(signif,"18e20 POT #nu","L");
202  leg->AddEntry(signif2,"6e20 POT #nu","L");
203  leg->AddEntry(signif3,"2e20 POT #nu","L");
204  leg->Draw("same");
205 
206  Simulation();
207 
208  gPad->Print("plots/signif.ps");
209  gPad->Print("plots/signif.C");
210 
211 }
void Simulation()
void SetDmsq21(const T &dmsq21) override
Definition: OscCalcPMNS.h:35
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void SetTh12(const T &th12) override
Definition: OscCalcPMNS.h:37
void SetTh23(const T &th23) override
Definition: OscCalcPMNS.h:39
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:18
Adapt the PMNS calculator to standard interface.
Definition: StanTypedefs.h:28
string fnameSwap
Definition: demo4.py:6
void SetL(double L) override
Definition: OscCalcPMNS.h:33
T sqrt(T number)
Definition: d0nt_math.hpp:156
const Cut kNumuFD
Definition: NumuCuts.h:53
osc::OscCalcDumb calc
void numu_sig_nonmax()
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const Binning kNumuEnergyBinning
Definition: Binnings.cxx:13
Spectrum Predict(osc::IOscCalc *calc) const override
Spectrum FakeData(double pot) const
Synonymous with AsimovData(). Retained for compatibility.
Definition: Spectrum.cxx:349
void SetdCP(const T &dCP) override
Definition: OscCalcPMNS.h:40
virtual std::unique_ptr< IFitSummary > Fit(osc::IOscCalcAdjustable *seed, SystShifts &bestSysts=junkShifts, const SeedList &seedPts=SeedList(), const std::vector< SystShifts > &systSeedPts={}, Verbosity verb=kVerbose) const
Master fitting method. Depends on FitHelper and FitHelperSeeded.
Definition: IFitter.cxx:69
void SetDmsq32(const T &dmsq32) override
Definition: OscCalcPMNS.h:36
const Var kCCE
Definition: NumuVars.h:21
#define pot
virtual void Go() override
Load all the registered spectra.
void SetTh13(const T &th13) override
Definition: OscCalcPMNS.h:38
loader
Definition: demo0.py:10
TLatex * prelim
Definition: Xsec_final.C:133
void SetRho(double rho) override
Definition: OscCalcPMNS.h:34
OStream cout
Definition: OStream.cxx:6
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
Prediction that just uses FD MC, with no extrapolation.
Float_t e
Definition: plot.C:35
loaderSwap
Definition: demo4.py:9
T asin(T number)
Definition: d0nt_math.hpp:60
Compare a single data spectrum to the MC + cosmics expectation.
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17
enum BeamMode string