saveContours_complete.C
Go to the documentation of this file.
1 #ifdef __CINT__
3  std::cout << "Sorry, you must run in compiled mode" << std::endl;
4 }
5 #else
6 
7 #include "../WorkshopIncludes.h"
8 
10 {
11  // Period 1
12  TFile* fPred_p1 = new TFile("../Predictions/Prediction_period1_systs_kCCE.root");
13  PredictionInterp* prediction_ugly_p1 = LoadFrom<PredictionInterp>(fPred_p1->GetDirectory("pred_ugly")).release();
14  PredictionInterp* prediction_bad_p1 = LoadFrom<PredictionInterp>(fPred_p1->GetDirectory("pred_bad")).release();
15  PredictionInterp* prediction_good_p1 = LoadFrom<PredictionInterp>(fPred_p1->GetDirectory("pred_good")).release();
16 
17  PredictionInterp* prediction_std_p1 = LoadFrom<PredictionInterp>(fPred_p1->GetDirectory("pred_std")).release();
18 
19  fPred_p1->Close();
20  // Period 2
21  TFile* fPred_p2 = new TFile("../Predictions/Prediction_period2_systs_kCCE.root");
22  PredictionInterp* prediction_ugly_p2 = LoadFrom<PredictionInterp>(fPred_p2->GetDirectory("pred_ugly")).release();
23  PredictionInterp* prediction_bad_p2 = LoadFrom<PredictionInterp>(fPred_p2->GetDirectory("pred_bad")).release();
24  PredictionInterp* prediction_good_p2 = LoadFrom<PredictionInterp>(fPred_p2->GetDirectory("pred_good")).release();
25 
26  PredictionInterp* prediction_std_p2 = LoadFrom<PredictionInterp>(fPred_p2->GetDirectory("pred_std")).release();
27 
28  fPred_p2->Close();
29  // Epoch 3b
30  TFile* fPred_3b = new TFile("../Predictions/Prediction_epoch3b_systs_kCCE.root");
31  PredictionInterp* prediction_ugly_3b = LoadFrom<PredictionInterp>(fPred_3b->GetDirectory("pred_ugly")).release();
32  PredictionInterp* prediction_bad_3b = LoadFrom<PredictionInterp>(fPred_3b->GetDirectory("pred_bad")).release();
33  PredictionInterp* prediction_good_3b = LoadFrom<PredictionInterp>(fPred_3b->GetDirectory("pred_good")).release();
34 
35  PredictionInterp* prediction_std_3b = LoadFrom<PredictionInterp>(fPred_3b->GetDirectory("pred_std")).release();
36 
37  fPred_3b->Close();
38  // Epoch 3c
39  TFile* fPred_3c = new TFile("../Predictions/Prediction_period1_systs_kCCE.root");
40  PredictionInterp* prediction_ugly_3c = LoadFrom<PredictionInterp>(fPred_3c->GetDirectory("pred_ugly")).release();
41  PredictionInterp* prediction_bad_3c = LoadFrom<PredictionInterp>(fPred_3c->GetDirectory("pred_bad")).release();
42  PredictionInterp* prediction_good_3c = LoadFrom<PredictionInterp>(fPred_3c->GetDirectory("pred_good")).release();
43 
44  PredictionInterp* prediction_std_3c = LoadFrom<PredictionInterp>(fPred_3c->GetDirectory("pred_std")).release();
45 
46  fPred_3c->Close();
47 
48  // Now let's combine the different periods
49 
50  PredictionCombinePeriods* prediction_ugly = new PredictionCombinePeriods({
51  {prediction_ugly_p1, kSecondAnaPeriod1POT},
52  {prediction_ugly_p2, kSecondAnaPeriod2POT},
53  {prediction_ugly_3b, kSecondAnaEpoch3bPOT},
54  {prediction_ugly_3c, kSecondAnaEpoch3cPOT + kSecondAnaEpoch3dPOT} });
55 
56  PredictionCombinePeriods* prediction_bad = new PredictionCombinePeriods({
57  {prediction_bad_p1, kSecondAnaPeriod1POT},
58  {prediction_bad_p2, kSecondAnaPeriod2POT},
59  {prediction_bad_3b, kSecondAnaEpoch3bPOT},
60  {prediction_bad_3c, kSecondAnaEpoch3cPOT + kSecondAnaEpoch3dPOT} });
61 
62  PredictionCombinePeriods* prediction_good = new PredictionCombinePeriods({
63  {prediction_good_p1, kSecondAnaPeriod1POT},
64  {prediction_good_p2, kSecondAnaPeriod2POT},
65  {prediction_good_3b, kSecondAnaEpoch3bPOT},
66  {prediction_good_3c, kSecondAnaEpoch3cPOT + kSecondAnaEpoch3dPOT} });
67 
68  PredictionCombinePeriods* prediction_std = new PredictionCombinePeriods({
69  {prediction_std_p1, kSecondAnaPeriod1POT},
70  {prediction_std_p2, kSecondAnaPeriod2POT},
71  {prediction_std_3b, kSecondAnaEpoch3bPOT},
72  {prediction_std_3c, kSecondAnaEpoch3cPOT + kSecondAnaEpoch3dPOT} });
73 
74  // And we're ready to go
75 
76  osc::OscCalculatorPMNSOpt calc;
77  ana::ResetOscCalcToDefault(&calc); // This seeds to default (in particular max-mixing)
78 
79  calc.SetTh23(TMath::ASin(sqrt(0.403))); // Second analysis best fit
80  calc.SetDmsq32(2.67e-3); // Same
81 
82  std::vector<const ISyst*> systs = getAllNumuSysts();
83 
84  std::vector<const IFitVar*> oscpar = getNumuMarginalisedOscParam();
85  std::vector<const IFitVar*> margPars = getNumuMarginalisedOscParam();
86  margPars.push_back(&kFitDmSq32Scaled);
87 
88  Spectrum fakePred = prediction_std->Predict(&calc);
89  Spectrum fakeData = fakePred.FakeData(kSecondAnaPOT);
90 
91  // Split experiments data
92  Spectrum fakePred_ugly = prediction_ugly->Predict(&calc);
93  Spectrum fakeData_ugly = fakePred_ugly.FakeData(kSecondAnaPOT);
94 
95  Spectrum fakePred_bad = prediction_bad->Predict(&calc);
96  Spectrum fakeData_bad = fakePred_bad.FakeData(kSecondAnaPOT);
97 
98  Spectrum fakePred_good = prediction_good->Predict(&calc);
99  Spectrum fakeData_good = fakePred_good.FakeData(kSecondAnaPOT);
100 
101  // Experiments with constraints
102 
103  SingleSampleExperiment fakeExpt(prediction_std, fakeData);
104  std::vector<const IExperiment*> vConstraints;
105  vConstraints.push_back(&fakeExpt);
106  vConstraints.push_back(new SolarConstraints());
107  vConstraints.push_back(WorldReactorConstraint2016());
108  MultiExperiment ExptWconstr(vConstraints);
109 
110 
111  // Multiple exprts
112  SingleSampleExperiment fakeExpt_ugly(prediction_ugly, fakeData_ugly);
113  SingleSampleExperiment fakeExpt_bad(prediction_bad, fakeData_bad);
114  SingleSampleExperiment fakeExpt_good(prediction_good, fakeData_good);
115  /// Combine experiments
116  std::vector<const IExperiment*> vExpts;
117  vExpts.push_back(&fakeExpt_ugly);
118  vExpts.push_back(&fakeExpt_bad);
119  vExpts.push_back(&fakeExpt_good);
120  vExpts.push_back(new SolarConstraints());
121  vExpts.push_back(WorldReactorConstraint2016());
122 
123  MultiExperiment multiExpt(vExpts);
124 
125  // 2D Fits
126 
127  Surface* surface = new Surface( &ExptWconstr, &calc,
128  &kFitSinSqTheta23, 40, 0.3, 0.7,
129  &kFitDmSq32Scaled, 35, 2.3, 3.0, oscpar, systs);
130 
131  Surface* surfaceExtra = new Surface( &multiExpt, &calc,
132  &kFitSinSqTheta23, 40, 0.3, 0.7,
133  &kFitDmSq32Scaled, 35, 2.3, 3.0, oscpar, systs);
134 
135  // Plot
136  TCanvas* c = new TCanvas();
137  c->SetLeftMargin(0.12);
138  c->SetBottomMargin(0.15);
139 
140  surface->DrawBestFit(kBlue);
141  surface->DrawContour( Gaussian90Percent2D(*surface), kSolid, kBlue);
142 
143  surfaceExtra->DrawBestFit(kRed);
144  surfaceExtra->DrawContour( Gaussian90Percent2D(*surfaceExtra), kDashed, kRed);
145 
146  c->SaveAs("../Contours/Complete.root");
147  c->SaveAs("../Contours/Complete.pdf");
148 
149  c->Clear();
150 
151  // Get 1D marginalised plots
152  TH1* profileSinsq = Profile(&ExptWconstr, &calc,
153  &kFitSinSqTheta23, 40, 0.3, 0.7, -1,
154  margPars, systs );
155 
156  TH1* profileSinsqMulti = Profile(&multiExpt, &calc,
157  &kFitSinSqTheta23, 40, 0.3, 0.7, -1,
158  margPars, systs );
159 
160  profileSinsq->SetLineColor(kBlue);
161  profileSinsq->Draw("hist");
162  profileSinsqMulti->SetLineColor(kRed);
163  profileSinsqMulti->SetLineStyle(7);
164  profileSinsqMulti->Draw("hist same");
165 
166  std::cout << "Max-mix rej 1 : " << sqrt(profileSinsq->Interpolate(0.5)) << std::endl;
167  std::cout << "Max-mix rej Multi: " << sqrt(profileSinsqMulti->Interpolate(0.5)) << std::endl;
168 
169  c->SaveAs("../Contours/1D_Complete.root");
170  c->SaveAs("../Contours/1D_Complete.pdf");
171  c->Close();
172  c->Close();
173 
174 } // End of function
175 
176 
177 #endif
Implements systematic errors by interpolation between shifted templates.
enum BeamMode kRed
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
const double kSecondAnaPeriod2POT
Definition: Exposures.h:74
Constraints on the parameters and from solar experiments.
T sqrt(T number)
Definition: d0nt_math.hpp:156
const double kSecondAnaEpoch3dPOT
Definition: Exposures.h:77
const double kSecondAnaEpoch3bPOT
Definition: Exposures.h:75
TGraph * Profile(const IExperiment *expt, osc::IOscCalcAdjustable *calc, const IFitVar *v, int nbinsx, double minx, double maxx, double input_minchi, const std::vector< const IFitVar * > &profVars, const std::vector< const ISyst * > &profSysts, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts, std::map< const IFitVar *, TGraph * > &profVarsMap, std::map< const ISyst *, TGraph * > &profSystsMap, MinuitFitter::FitOpts opts)
scan in one variable, profiling over all others
Definition: Fit.cxx:48
const double kSecondAnaPeriod1POT
Definition: Exposures.h:73
osc::OscCalcDumb calc
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
void saveContours_complete()
Spectrum FakeData(double pot) const
Synonymous with AsimovData(). Retained for compatibility.
Definition: Spectrum.cxx:349
const ReactorExperiment * WorldReactorConstraint2016()
Updated value for SecondAna based on the latest PDG.
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
OStream cout
Definition: OStream.cxx:6
Combine multiple component experiments.
TH2 * Gaussian90Percent2D(const FrequentistSurface &s)
Up-value surface for 90% confidence in 2D in gaussian approximation.
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
const double kSecondAnaEpoch3cPOT
Definition: Exposures.h:76
Sum MC predictions from different periods scaled according to data POT targets.
Float_t e
Definition: plot.C:35
const double kSecondAnaPOT
Definition: Exposures.h:87
enum BeamMode kBlue
Spectrum Predict(osc::IOscCalc *calc) const override
Compare a single data spectrum to the MC + cosmics expectation.
std::vector< const IFitVar * > getNumuMarginalisedOscParam()
Definition: FitVarsNumu.cxx:9