makeMatrixElementSurface.C
Go to the documentation of this file.
1 #include <cmath>
2 #include <iostream>
3 #include <iomanip>
4 #include <string>
5 #include <vector>
6 
10 #include "CAFAna/Analysis/Plots.h"
11 #include "CAFAna/Analysis/Style.h"
19 #include "NuXAna/Systs/NusSysts.h"
20 #include "CAFAna/Systs/Systs.h"
21 #include "CAFAna/Vars/FitVars.h"
25 #include "OscLib/OscCalcSterile.h"
26 #include "CAFAna/Core/Progress.h"
27 #include "CAFAna/Core/Utilities.h"
28 
29 #include "TCanvas.h"
30 #include "TFile.h"
31 #include "TH1.h"
32 #include "TH2.h"
33 #include "TMath.h"
34 #include "TGaxis.h"
35 #include "TLegend.h"
36 #include "TLatex.h"
37 #include "TStyle.h"
38 #include "TPaveText.h"
39 #include "TLegendEntry.h"
40 #include "TMarker.h"
41 
42 
43 using namespace ana;
44 
45 
46 void TexInfo();
47 
48 TLegend* MakeLegend(double x1, double y1, double x2, double y2, double textsize = 0.03)
49 {
50  TLegend* leg = new TLegend(x1, y1, x2, y2);
51  leg->SetBorderSize(0);
52  leg->SetFillColor(kWhite);
53  leg->SetFillStyle(0);
54  leg->SetFillStyle(0);
55  leg->SetTextFont(42);
56  leg->SetTextSize(textsize);
57  return leg;
58 }
59 
60 
61 void SigmaLines(double xmin, double xmax);
62 
63 void PlotText(const double xpos, const double start, const double step);
64 
65 // Reset calculator values after each fit to preset defaults
67 
69 {
70  TH1::AddDirectory(0);
71 
72 
73  TFile input(inputFilename.c_str(), "READ"); // Load in MC prediction
74 
75  // Get the nus systematics
76  std::vector<const ISyst*> systs = getAllNusSysts();
77 
78  // Load in the NuMI data spectrum
79  Spectrum* intime = LoadFromFile<Spectrum>("$NOVA_ANA/steriles/Ana01/fout_nus_ana01_box_opening_restricted.root","spec_nus_cale_numi").release();
80 
81  TH1 *hSpectrum = intime->ToTH1(intime->POT());
82  Spectrum data_spec(hSpectrum, intime->POT(), intime->Livetime());
83  std::cout << "Data spectrum: " << data_spec.POT()
84  <<" POT and "
85  << data_spec.ToTH1(intime->POT())->GetSumOfWeights()
86  << " events" << std::endl;
87 
88  // Load in the out-of-NuMI-time weighted spectrum (cosmics)
89  Spectrum* cosmic_scale = LoadFromFile<Spectrum>("$NOVA_ANA/steriles/Ana01/AnaResults.root","cosmic__CalE").release();
90  TH1D* hcosmic = (TH1D*)cosmic_scale->ToTH1(cosmic_scale->Livetime(), kLivetime);
91  Spectrum cosmic(hcosmic, intime->POT(), intime->Livetime());
92 
93  // Load the MC Prediction
94  std::unique_ptr<PredictionInterp> pred_p1p2e3b =
95  PredictionInterp::LoadFrom((TDirectory*)input.Get("nus_pred_p1p2e3b"));
96 
97  std::unique_ptr<PredictionInterp> pred_e3c =
98  PredictionInterp::LoadFrom((TDirectory*)input.Get("nus_pred_e3c"));
99 
100  const double pot_p1p2e3b =
104 
105  const double pot_e3ce3d = kSecondAnaEpoch3cPOT + kSecondAnaEpoch3dPOT;
106 
107  // Combine period1/period2/epoch3b + epoch3c, weighted appropriately
109  {&(*pred_p1p2e3b),pot_p1p2e3b},
110  {&(*pred_e3c),pot_e3ce3d}
111  });
112 
113  // Make a default sterile calculator and set values
115  calc4f->SetNFlavors(4);
116  ResetAngles(calc4f); // paranoia
117 
118  // Output file
119  TFile* rootOut = new TFile("matrixElements_osc_syst_incl_fix.root","recreate");
120 
121  // Take 68% band from PDG LIVE 2015: pdglive.lbl.gov
122  // #sin^{2}#theta_{23} = 0.514 + 0.055 - 0.056
123  // #asin(sqrt(0.514)) + #asin(sqrt(0.541 + 0.055)) - #asin(sqrt(0.541 + 0.056)) in degrees
124  // This is the Gaussian constraint on #theta_{23}
126  45.8,
127  3.2);
128 
129  CountingExperiment nus_expt(&pred, data_spec, cosmic);
130  MultiExperiment nus_multi({&th23Constraint, &nus_expt});
131 
132 
133  // Vector of parameters to fit: #theta_{23}, #theta_{34}, #theta_{24}
134  std::vector<const IFitVar*> fitvars;
135  fitvars.push_back(&kFitTheta34InDegreesSterile);
136  fitvars.push_back(&kFitTheta24InDegreesSterile);
137  fitvars.push_back(&kFitTheta23InDegreesSterile);
138  fitvars.push_back(&kFitDelta24InPiUnitsSterile);
139 
140 
141  const Color_t kStatsFitCol = kAzure+2;
142  const Color_t kSystsFitCol = kBlack;
143 
144  // Fit |Umu4|^2 vs |Utau4|^2
145  int nUmu4 = 200;
146  double minUmu4 = 0;
147  double maxUmu4 = 0.35;
148  int nUtau4 = 200;
149  double minUtau4 = 0;
150  double maxUtau4 = 0.30;
151 
152  Progress* prog = new Progress("hUtau4Umu4");
153  TH2F* hUmu4Utau4 = ExpandedHistogram(";|U_{#mu4}|^{2};|U_{#tau4}|^{2}", nUmu4, minUmu4, maxUmu4, false, nUtau4, minUtau4, maxUtau4, false);
154  for (int iUmu4 = 0; iUmu4 < nUmu4; ++iUmu4)
155  {
156  for (int iUtau4 = 0; iUtau4 < nUtau4; ++iUtau4)
157  {
158  double Umu4 = hUmu4Utau4->GetXaxis()->GetBinCenter(iUmu4+1);
159  double Utau4 = hUmu4Utau4->GetYaxis()->GetBinCenter(iUtau4+1);
160  double th24 = asin( sqrt(Umu4) );
161  double th34 = asin( sqrt(Utau4/(1 - Umu4)) );
162  ResetAngles(calc4f);
163  kFitTheta24Sterile.SetValue(calc4f, th24);
164  kFitTheta34Sterile.SetValue(calc4f, th34);
166  SystShifts systSeed = SystShifts::Nominal();
167  double chi2 = fitter.Fit(calc4f, systSeed, IFitter::kQuiet)->EvalMetricVal();
168  hUmu4Utau4->SetBinContent(iUmu4+1, iUtau4+1, chi2);
169  prog->SetProgress( (iUmu4*nUmu4 + iUtau4)/double(nUtau4*nUmu4) );
170  }
171  }
172 
173  double minchi2 = 1e10;
174  for (int iUtau4 = 0; iUtau4 < nUtau4; ++iUtau4)
175  {
176  for (int iUmu4 = 0; iUmu4 < nUmu4; ++iUmu4)
177  {
178  double chi2 = hUmu4Utau4->GetBinContent(iUtau4+1, iUmu4+1);
179  if (chi2 < minchi2) minchi2 = chi2;
180  }
181  }
182  for (int iUtau4 = 0; iUtau4 < nUtau4; ++iUtau4)
183  {
184  for (int iUmu4 = 0; iUmu4 < nUmu4; ++iUmu4)
185  {
186  double chi2 = hUmu4Utau4->GetBinContent(iUtau4+1, iUmu4+1);
187  chi2 -= minchi2;
188  hUmu4Utau4->SetBinContent(iUtau4+1, iUmu4+1, chi2);
189  }
190  }
191  rootOut->WriteTObject(hUmu4Utau4, "hUmu4Utau4");
192 
193  rootOut->Close();
194  std::cout << "All done!" << std::endl;
195 }
196 
197 
198 
199 
200 void PlotText(const double xpos, const double start, const double step)
201 {
202  TLatex* tex = new TLatex();
203  tex->SetNDC();
204  tex->SetTextFont(42);
205  tex->SetTextSize(0.037);
206  tex->SetTextAlign(11);
207  tex->DrawLatex(xpos, start, "NOvA 6.05#times10^{20} POT-equiv.");
208  tex->DrawLatex(xpos, start - 1*step, "#Deltam^{2}_{32} = 2.44#times10^{-3} eV^{2}");
209  tex->DrawLatex(xpos, start - 2*step, "#theta_{13} = 8.5#circ, ^{}#theta_{23} = 45#circ,");
210  tex->DrawLatex(xpos, start - 3*step, "#Deltam^{2}_{41} = 0.5 eV^{2}");
211 }
212 
213 
214 void TexInfo()
215 {
216  const double xtex = 0.76;
217  const double y68 = 0.28;
218  const double y90 = 0.55;
219 
220  TLatex *tex = new TLatex();
221  tex->SetNDC();
222  tex->SetTextFont(42);
223  tex->SetTextSize(0.037);
224  tex->SetTextAlign(11);
225  tex->DrawLatex(xtex, y68, "68% C.L."); // 1-sigma
226  tex->DrawLatex(xtex, y90, "90% C.L."); // 2-sigma
227 }
228 
229 
230 void SigmaLines(double xmin, double xmax)
231 {
232  TLine* line68 = new TLine();
233  TLine* line90 = new TLine();
234  // Fixes horizontal effects at max canvas value
235  TLine* line05 = new TLine();
236  // Dashed line for 1 and 2 sigma lines
237  line68->SetLineStyle(kDashed);
238  line90->SetLineStyle(kDashed);
239 
240  line68->SetLineWidth(2);
241  line90->SetLineWidth(2);
242  line05->SetLineWidth(2);
243  line68->DrawLine(xmin, 1, xmax, 1);
244  line90->DrawLine(xmin, 2.71, xmax, 2.71);
245  line05->DrawLine(xmin, 5, xmax, 5);
246 }
247 
248 
250 {
252  calc->SetDm(4, 0.5);
253  calc->SetAngle(2, 3, M_PI/4); // 45 degress, maximal mixing
254  calc->SetAngle(3, 4, 0.); // 0 degrees
255  calc->SetAngle(2, 4, 0.); // 0 degrees
256  return;
257 }
TCut intime("tave>=217.0 && tave<=227.8")
void makeMatrixElementSurface(std::string inputFilename)
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void ResetAngles(osc::OscCalcSterile *calc)
std::map< std::string, double > xmax
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
const double kSecondAnaPeriod2POT
Definition: Exposures.h:74
Float_t y1[n_points_granero]
Definition: compare.C:5
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
A simple Gaussian constraint on an arbitrary IFitVar.
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
Float_t x1[n_points_granero]
Definition: compare.C:5
T sqrt(T number)
Definition: d0nt_math.hpp:156
const double kSecondAnaEpoch3dPOT
Definition: Exposures.h:77
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
const double kSecondAnaEpoch3bPOT
Definition: Exposures.h:75
Adapt the PMNS_Sterile calculator to standard interface.
const FitTheta23InDegreesSterile kFitTheta23InDegreesSterile
const FitTheta24Sterile kFitTheta24Sterile
const FitTheta34Sterile kFitTheta34Sterile
const double kSecondAnaPeriod1POT
Definition: Exposures.h:73
const FitTheta24InDegreesSterile kFitTheta24InDegreesSterile
static SystShifts Nominal()
Definition: SystShifts.h:34
osc::OscCalcDumb calc
double chi2()
#define M_PI
Definition: SbMath.h:34
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
TH2F * ExpandedHistogram(const std::string &title, int nbinsx, double xmin, double xmax, bool xlog, int nbinsy, double ymin, double ymax, bool ylog)
Internal helper for Surface and FCSurface.
Definition: Utilities.cxx:149
const FitDelta24InPiUnitsSterile kFitDelta24InPiUnitsSterile
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
std::vector< const ISyst * > getAllNusSysts()
Get a vector of all the nus group systs.
Definition: NusSysts.cxx:202
Sum up livetimes from individual cosmic triggers.
void TexInfo()
void ResetSterileCalcToDefault(osc::OscCalcSterile *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:76
double POT() const
Definition: Spectrum.h:227
TLegend * MakeLegend(double x1, double y1, double x2, double y2, double textsize=0.03)
OStream cout
Definition: OStream.cxx:6
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
void PlotText(const double xpos, const double start, const double step)
Combine multiple component experiments.
void SetAngle(int i, int j, double th)
static std::unique_ptr< PredictionInterp > LoadFrom(TDirectory *dir, const std::string &name)
void SetProgress(double frac)
Update the progress fraction between zero and one.
Definition: Progress.cxx:41
void SetDm(int i, double dm)
TLatex * tex
Definition: f2_nu.C:499
A simple ascii-art progress bar.
Definition: Progress.h:9
const FitTheta34InDegreesSterile kFitTheta34InDegreesSterile
const double kSecondAnaEpoch3cPOT
Definition: Exposures.h:76
Sum MC predictions from different periods scaled according to data POT targets.
void SigmaLines(double xmin, double xmax)
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:230
Compare a data spectrum to MC expectation using only the event count.
T asin(T number)
Definition: d0nt_math.hpp:60
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17
enum BeamMode string