demo_CPT.C
Go to the documentation of this file.
1 #include "CAFAna/Vars/Vars.h"
10 #include "OscLib/OscCalcCPT.h"
11 #include "TCanvas.h"
12 #include "TGaxis.h"
13 #include "TMath.h"
14 
15 using namespace ana;
16 
17 void demo_CPT()
18 {
19 
20  const std::string fname =
21  "$NOVA_PROD/mc/FA14-10-28/genie/fd/caf/010000/1000000/*_fhc_nonswap_*";
22  const std::string fnameSwap =
23  "$NOVA_PROD/mc/FA14-10-28/genie/fd/caf/010000/1000000/*_fhc_swap_*";
24  const std::string fnameRHC =
25  "$NOVA_PROD/mc/FA14-10-28/genie/fd/caf/010000/1000000/*_rhc_nonswap_*";
26  const std::string fnameRHCSwap =
27  "$NOVA_PROD/mc/FA14-10-28/genie/fd/caf/010000/1000000/*_rhc_swap_*";
28 
29  SpectrumLoader loader(fname);
30  SpectrumLoader loaderSwap(fnameSwap);
31  SpectrumLoader loaderRHC(fnameRHC);
32  SpectrumLoader loaderRHCSwap(fnameRHCSwap);
33 
34  PredictionNoExtrap predCC( loader, loaderSwap,
35  "Simple Energy (GeV)", kNumuEnergyBinning,
36  kCCE, kNumuNonQEFD );
37  PredictionNoExtrap predQE( loader, loaderSwap,
38  "QE Energy (GeV)", kNumuEnergyBinning,
39  kQEE, kNumuQEFD );
40  PredictionNoExtrap predCCRHC(loaderRHC, loaderRHCSwap,
41  "Simple Energy (GeV)", kNumuEnergyBinning,
42  kCCE, kNumuNonQEFD );
43  PredictionNoExtrap predQERHC(loaderRHC, loaderRHCSwap,
44  "QE Energy (GeV)", kNumuEnergyBinning,
45  kQEE, kNumuQEFD );
46 
47 
48  loader.Go();
49  loaderSwap.Go();
50  loaderRHC.Go();
51  loaderRHCSwap.Go();
52 
53 
54  const double pot = 18e20;
55  const double potrhc = 18e20;
56 
57  //Set up the CPT-violating oscillation calculators
58  osc::OscCalcCPT calcR;
59  calcR.SetL(810);
60  calcR.SetRho(2.75);
61  calcR.SetDmsq21(7.59e-5);
62  calcR.SetDmsq32(2.41e-3);
63  calcR.SetTh12(.601); //sin2(2t12) = .870
64  calcR.SetTh13(.1567); //sin2(2t13) = .095
65  calcR.SetTh23(TMath::Pi()*.2141); //sin2(2t23) = .95
66 
67  osc::OscCalcCPT calcB;
68  calcB.SetL(810);
69  calcB.SetRho(2.75);
70  calcB.SetDmsq21(7.59e-5);
71  calcB.SetDmsq32(2.41e-3, osc::ENuSign::kNu);
72  calcB.SetDmsq32(3.00e-3, osc::ENuSign::kNuBar);
73  calcB.SetTh12(.601); //sin2(2t12) = .870
74  calcB.SetTh13(.1567); //sin2(2t13) = .095
75  calcB.SetTh23(TMath::Pi()*.2141, osc::ENuSign::kNu); //sin2(2t23) = .95
76  calcB.SetTh23(TMath::Pi()*.1867, osc::ENuSign::kNuBar); //sin2(2t23) = .85
77 
78  //Make fake data
79  Spectrum obsCCR = predCC.Predict(&calcR);
80  Spectrum obsQER = predQE.Predict(&calcR);
81  Spectrum obsCCRHCR = predCCRHC.Predict(&calcR);
82  Spectrum obsQERHCR = predQERHC.Predict(&calcR);
83 
84  Spectrum fakeCCR = obsCCR.FakeData(pot);
85  Spectrum fakeQER = obsQER.FakeData(pot);
86  Spectrum fakeCCRHCR = obsCCRHCR.FakeData(potrhc);
87  Spectrum fakeQERHCR = obsQERHCR.FakeData(potrhc);
88 
89 
90  Spectrum obsCCB = predCC.Predict(&calcB);
91  Spectrum obsQEB = predQE.Predict(&calcB);
92  Spectrum obsCCRHCB = predCCRHC.Predict(&calcB);
93  Spectrum obsQERHCB = predQERHC.Predict(&calcB);
94 
95  Spectrum fakeCCB = obsCCB.FakeData(pot);
96  Spectrum fakeQEB = obsQEB.FakeData(pot);
97  Spectrum fakeCCRHCB = obsCCRHCB.FakeData(potrhc);
98  Spectrum fakeQERHCB = obsQERHCB.FakeData(potrhc);
99 
100  //Make two experiments
101  SingleSampleExperiment exptCCR(&predCC, fakeCCR);
102  SingleSampleExperiment exptQER(&predQE, fakeQER);
103  SingleSampleExperiment exptCCRHCR(&predCCRHC, fakeCCRHCR);
104  SingleSampleExperiment exptQERHCR(&predQERHC, fakeQERHCR);
105 
106  std::vector<const IExperiment*> exptsR;
107 
108  exptsR.push_back(&exptCCR);
109  exptsR.push_back(&exptQER);
110  exptsR.push_back(&exptCCRHCR);
111  exptsR.push_back(&exptQERHCR);
112 
113  MultiExperiment exptR(exptsR);
114 
115 
116  SingleSampleExperiment exptCCB(&predCC, fakeCCB);
117  SingleSampleExperiment exptQEB(&predQE, fakeQEB);
118  SingleSampleExperiment exptCCRHCB(&predCCRHC, fakeCCRHCB);
119  SingleSampleExperiment exptQERHCB(&predQERHC, fakeQERHCB);
120 
121  std::vector<const IExperiment*> exptsB;
122 
123  exptsB.push_back(&exptCCB);
124  exptsB.push_back(&exptQEB);
125  exptsB.push_back(&exptCCRHCB);
126  exptsB.push_back(&exptQERHCB);
127 
128  MultiExperiment exptB(exptsB);
129 
130  //Draw the surfaces
131  TGaxis::SetMaxDigits(3);
132 
133  const int res = 15; //10 is fast but ugly, 100 is slow and pretty
134  //pick something in between
135  const double xmin = -.15;
136  const double xmax = .25;
137  const double ymin = -1.5e-3;
138  const double ymax = 1e-3;
139 
140  TCanvas* c1 = new TCanvas("c1", "cptPlot", 1024, 768);
141 
142 
145  FitDeltaCPT kFitDeltaCPTSinSq2Theta23B(&sdtB);
146  FitSigmaCPT kFitSigmaCPTSinSq2Theta23B(&sdtB);
147  FitDeltaCPT kFitDeltaCPTDmSq32NormB(&sdmB);
148  FitSigmaCPT kFitSigmaCPTDmSq32NormB(&sdmB);
149 
150  FrequentistSurface surfB2(&exptB, &calcB,
151  &kFitDeltaCPTSinSq2Theta23B, res, xmin, xmax,
152  &kFitDeltaCPTDmSq32NormB, res, ymin, ymax,
153  {&kFitSigmaCPTSinSq2Theta23B,&kFitSigmaCPTDmSq32NormB}
154  );
155 
156  surfB2.DrawBestFit(kBlue);
157  surfB2.DrawContour(Gaussian68Percent2D(surfB2), 1, kBlue);
158  surfB2.DrawContour(Gaussian2Sigma2D(surfB2), 2, kBlue);
159  surfB2.DrawContour(Gaussian3Sigma2D(surfB2), 3, kBlue);
160 
161 
162  TLine *xaxis = new TLine(xmin,0,xmax,0);
163  xaxis->SetLineColor(kGray);
164  xaxis->SetLineStyle(7);
165  xaxis->Draw();
166 
167  TLine *yaxis = new TLine(0,ymin,0,ymax);
168  yaxis->SetLineColor(kGray);
169  yaxis->SetLineStyle(7);
170  yaxis->Draw();
171 
172 
175  FitDeltaCPT kFitDeltaCPTSinSq2Theta23R(&sdtR);
176  FitSigmaCPT kFitSigmaCPTSinSq2Theta23R(&sdtR);
177  FitDeltaCPT kFitDeltaCPTDmSq32NormR(&sdmR);
178  FitSigmaCPT kFitSigmaCPTDmSq32NormR(&sdmR);
179 
180  FrequentistSurface surfR2(&exptR, &calcR,
181  &kFitDeltaCPTSinSq2Theta23R, res, xmin, xmax,
182  &kFitDeltaCPTDmSq32NormR, res, ymin, ymax,
183  {&kFitSigmaCPTSinSq2Theta23R,&kFitSigmaCPTDmSq32NormR}
184  );
185 
186  surfR2.DrawBestFit(kRed);
187  surfR2.DrawContour(Gaussian68Percent2D(surfR2), 1, kRed);
188  surfR2.DrawContour(Gaussian2Sigma2D(surfR2), 2, kRed);
189  surfR2.DrawContour(Gaussian3Sigma2D(surfR2), 3, kRed);
190 
191 
192  c1->SaveAs("CPT-sesitivity.pdf");
193 
194 }
void SetDmsq21(const double &dmsq21) override
Definition: OscCalcCPT.h:66
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::map< std::string, double > xmax
void SetDmsq32(const double &dmsq32) override
Definition: OscCalcCPT.h:68
string fnameSwap
Definition: demo4.py:6
TH2 * Gaussian68Percent2D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 2D in gaussian approximation.
std::string yaxis
TH2 * Gaussian2Sigma2D(const FrequentistSurface &s)
Up-value surface for 2 sigma confidence in 2D in gaussian approximation.
std::string xaxis
const FitSinSq2Theta23CPT kFitSinSq2Theta23bar(osc::ENuSign::kNuBar)
Definition: FitVarsCPT.h:137
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
void SetTh12(const double &th12) override
Definition: OscCalcCPT.h:70
Double_t ymax
Definition: plot.C:25
Log-likelihood scan across two parameters.
void SetTh13(const double &th13) override
Definition: OscCalcCPT.h:72
void DrawBestFit(Color_t color, Int_t marker=kFullCircle) const
Draw the best fit point.
Definition: ISurface.cxx:32
void SetRho(double rho) override
Definition: OscCalcCPT.h:64
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
const FitDmSq32CPTHierarchy kFitDmSq32nuNorm(osc::ENuSign::kNu, false)
Definition: FitVarsCPT.h:89
const Var kCCE
Definition: NumuVars.h:21
#define pot
TH2 * Gaussian3Sigma2D(const FrequentistSurface &s)
Up-value surface for 3 sigma confidence in 2D in gaussian approximation.
const FitSinSq2Theta23CPT kFitSinSq2Theta23nu(osc::ENuSign::kNu)
Definition: FitVarsCPT.h:136
virtual void Go() override
Load all the registered spectra.
const Cut kNumuNonQEFD
Definition: NumuCuts.h:67
loader
Definition: demo0.py:10
void SetTh23(const double &th23) override
Definition: OscCalcCPT.h:74
void demo_CPT()
Definition: demo_CPT.C:17
Combine multiple component experiments.
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
void SetL(double L) override
Definition: OscCalcCPT.h:62
const Var kQEE
Energy estimator for quasielastic CC events.
Definition: NumuVars.cxx:28
c1
Definition: demo5.py:24
Double_t ymin
Definition: plot.C:24
const Cut kNumuQEFD
Definition: NumuCuts.h:65
Prediction that just uses FD MC, with no extrapolation.
Float_t e
Definition: plot.C:35
loaderSwap
Definition: demo4.py:9
enum BeamMode kBlue
const FitDmSq32CPTHierarchy kFitDmSq32barNorm(osc::ENuSign::kNuBar, false)
Definition: FitVarsCPT.h:90
Compare a single data spectrum to the MC + cosmics expectation.
enum BeamMode string