saveContours_systs.C
Go to the documentation of this file.
1 #ifdef __CINT__
2 void saveContours_systs(){
3  std::cout << "Sorry, you must run in compiled mode" << std::endl;
4 }
5 #else
6 
7 #include "../WorkshopIncludes.h"
8 
10 {
11  TFile* fPred_addExpt = new TFile("../Predictions/Prediction_epoch3b_someSysts_kCCE.root");
12  PredictionInterp* prediction_ugly = LoadFrom<PredictionInterp>(fPred_addExpt->GetDirectory("pred_ugly")).release();
13  PredictionInterp* prediction_bad = LoadFrom<PredictionInterp>(fPred_addExpt->GetDirectory("pred_bad")).release();
14  PredictionInterp* prediction_good = LoadFrom<PredictionInterp>(fPred_addExpt->GetDirectory("pred_good")).release();
15 
16  PredictionInterp* prediction_std = LoadFrom<PredictionInterp>(fPred_addExpt->GetDirectory("pred_std")).release();
17 
18  fPred_addExpt->Close();
19 
20  osc::OscCalculatorPMNSOpt calc;
21  ana::ResetOscCalcToDefault(&calc); // This seeds to default (in particular max-mixing)
22 
23  calc.SetTh23(TMath::ASin(sqrt(0.403))); // Second analysis best fit
24  calc.SetDmsq32(2.67e-3); // Same
25 
26  //std::vector<const ISyst*> systs = getAllNumuSysts(); -> contains all the numu systematics
27  std::vector<const ISyst*> someSysts = {&kNumuNCScaleSyst, &kSAMuEnergyScaleSyst,
28  &kFDSAMuEnergyScaleSyst, &kSARelHadESyst }; // short subset to run faster
29 
30  // std::vector<const IFitVar*> oscpar = getNumuMarginalisedOscParam(); --> This gives all the osc param
31  std::vector<const IFitVar*> someOscpar = {&kFitSinSq2Theta13}; // to seeep up, th13
32 
33  Spectrum fakePred = prediction_std->Predict(&calc);
34  Spectrum fakeData = fakePred.FakeData(kSecondAnaPOT);
35 
36  // Split experiments data
37  Spectrum fakePred_ugly = prediction_ugly->Predict(&calc);
38  Spectrum fakeData_ugly = fakePred_ugly.FakeData(kSecondAnaPOT);
39 
40  Spectrum fakePred_bad = prediction_bad->Predict(&calc);
41  Spectrum fakeData_bad = fakePred_bad.FakeData(kSecondAnaPOT);
42 
43  Spectrum fakePred_good = prediction_good->Predict(&calc);
44  Spectrum fakeData_good = fakePred_good.FakeData(kSecondAnaPOT);
45 
46  // Experiments with constraints
47 
48  SingleSampleExperiment fakeExpt(prediction_std, fakeData);
49  std::vector<const IExperiment*> vConstraints;
50  vConstraints.push_back(&fakeExpt);
51  vConstraints.push_back(new SolarConstraints());
52  vConstraints.push_back(WorldReactorConstraint2016());
53  MultiExperiment ExptWconstr(vConstraints);
54 
55 
56  // Multiple exprts
57  SingleSampleExperiment fakeExpt_ugly(prediction_ugly, fakeData_ugly);
58  SingleSampleExperiment fakeExpt_bad(prediction_bad, fakeData_bad);
59  SingleSampleExperiment fakeExpt_good(prediction_good, fakeData_good);
60  /// Combine experiments
61  std::vector<const IExperiment*> vExpts;
62  vExpts.push_back(&fakeExpt_ugly);
63  vExpts.push_back(&fakeExpt_bad);
64  vExpts.push_back(&fakeExpt_good);
65  vExpts.push_back(new SolarConstraints());
66  vExpts.push_back(WorldReactorConstraint2016());
67 
68  MultiExperiment multiExpt(vExpts);
69 
70  // 2D Fits
71 
72  Surface* surface = new Surface( &ExptWconstr, &calc,
73  &kFitSinSqTheta23, 40, 0.3, 0.7,
74  &kFitDmSq32Scaled, 35, 2.3, 3.0, someOscpar, someSysts);
75 
76  Surface* surfaceExtra = new Surface( &multiExpt, &calc,
77  &kFitSinSqTheta23, 40, 0.3, 0.7,
78  &kFitDmSq32Scaled, 35, 2.3, 3.0, someOscpar, someSysts);
79 
80  // Plot
81  TCanvas* c = new TCanvas();
82  c->SetLeftMargin(0.12);
83  c->SetBottomMargin(0.15);
84 
85  surface->DrawBestFit(kBlue);
86  surface->DrawContour( Gaussian90Percent2D(*surface), kSolid, kBlue);
87 
88  surfaceExtra->DrawBestFit(kRed);
89  surfaceExtra->DrawContour( Gaussian90Percent2D(*surfaceExtra), kDashed, kRed);
90 
91  c->SaveAs("../Contours/ExtrapOscparAndSysts.root");
92  c->SaveAs("../Contours/ExtrapOscparAndSysts.pdf");
93 
94  c->Clear();
95 
96  // Get 1D marginalised plots
97  TH1* profileSinsq = Profile(&ExptWconstr, &calc,
98  &kFitSinSqTheta23, 40, 0.3, 0.7, -1,
99  {&kFitDmSq32Scaled, &kFitSinSq2Theta13}, someSysts );
100 
101  TH1* profileSinsqMulti = Profile(&multiExpt, &calc,
102  &kFitSinSqTheta23, 40, 0.3, 0.7, -1,
103  {&kFitDmSq32Scaled, &kFitSinSq2Theta13}, someSysts );
104 
105  profileSinsq->SetLineColor(kBlue);
106  profileSinsq->Draw("hist");
107  profileSinsqMulti->SetLineColor(kRed);
108  profileSinsqMulti->SetLineStyle(7);
109  profileSinsqMulti->Draw("hist same");
110 
111  std::cout << "Max-mix rej 1 : " << sqrt(profileSinsq->Interpolate(0.5)) << std::endl;
112  std::cout << "Max-mix rej Multi: " << sqrt(profileSinsqMulti->Interpolate(0.5)) << std::endl;
113 
114  c->SaveAs("../Contours/1D_SomeSysts.root");
115  c->SaveAs("../Contours/1D_SomeSysts.pdf");
116  c->Close();
117  c->Close();
118 
119 } // End of function
120 
121 
122 #endif
const NumuNCScaleSyst kNumuNCScaleSyst
Definition: NumuSysts.cxx:27
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
Constraints on the parameters and from solar experiments.
void saveContours_systs()
T sqrt(T number)
Definition: d0nt_math.hpp:156
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
osc::OscCalcDumb calc
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
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
Spectrum Predict(osc::IOscCalc *calc) const override
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 FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
Float_t e
Definition: plot.C:35
const double kSecondAnaPOT
Definition: Exposures.h:87
enum BeamMode kBlue
Compare a single data spectrum to the MC + cosmics expectation.