FCTutorial2020.C
Go to the documentation of this file.
1 // cafana core
2 #include "CAFAna/Core/Spectrum.h"
3 #include "CAFAna/Prediction/PredictionSyst3Flavor2020.h"
4 #include "CAFAna/Core/Progress.h"
5 
6 // inputs for fits
10 
11 // fitting machinery
12 #include "CAFAna/Fit/Fit.h"
13 #include "OscLib/IOscCalc.h"
14 #include "CAFAna/Analysis/Calcs.h"
15 #include "CAFAna/Vars/FitVars.h"
16 
17 // POT scaling constants
19 
20 // FC classes
21 #include "CAFAna/FC/FCCollection.h"
22 #include "CAFAna/FC/FCPoint.h"
23 
24 #include "TFile.h"
25 #include "TGraph.h"
26 #include "TH1D.h"
27 #include "TCanvas.h"
28 using namespace ana;
29 
30 void PlotPseudoexperiments(std::vector<TH1*> pse, std::string name);
31 void FCTutorial2020(int NPSE = 10, int nbins = 1, int jobid = 0, int procid = 0)
32 {
33  auto calc = DefaultOscCalc();
34 
35  // load predictions
36  bool isFakeData = true;
37  int nquantiles = 4;
38 
40  calc,
41  "fhc",
42  1,
43  {},
44  isFakeData);
45 
46 
47  // load gaussian profiles from FCInput
48  TFile * fcinput = TFile::Open("FCTutorial2020_FCInput.root");
49  TGraph * profile_dmsq = (TGraph*) fcinput->Get("NH_DmSq32" );
50  TGraph * profile_th13 = (TGraph*) fcinput->Get("NH_SinSq2Theta13");
51  fcinput->Close();
52 
53  // object for saving results
54  FCCollection fccol;
55 
56  // range in \sin^2\theta_{23}
57  double min_th23 = 0.3;
58  double max_th23 = 0.7;
59  std::vector<double> th23_bins;
60 
61  // given number of bins, populate values of \sin^2\theta_{23} to test
62  for(auto ith23 = 0; ith23 < nbins; ith23++)
63  th23_bins.push_back(min_th23 + (max_th23 - min_th23) * ith23 / nbins );
64 
65  // progress bar
66  Progress * prog = 0;
67 
68  // loop over all null hypotheses
69  for(auto ith23 = 0u; ith23 < th23_bins.size(); ith23++) {
70  double th23 = th23_bins[ith23];
71  // get parameters with which to "throw" experiments
72  // from gaussian profiles
73 
74  // null hypothesis also includes these parameters
75  // so that they are "profiled" out
76  double dmsq = profile_dmsq->Eval(th23);
77  double th13 = profile_th13->Eval(th23);
78 
79  // just take default value of deltacp
81 
82  // set calculator
83  kFitDmSq32 .SetValue(calc, dmsq );
85 
86  prog = new Progress(TString::Format("Throwing %d experiments at SinSqTheta23 = %.2f", NPSE, th23).Data());
87  // loop over number of pseudoexperiments
88  for(auto n = 0; n < NPSE; n++) {
89  prog->SetProgress(double(n+1) / NPSE);
90 
91  // make mock data
92  // adds poisson fluctuations to nominal prediction at calc and
93  // scales to data POT
94  Spectrum mock_data = pred_numu->Predict(calc).MockData(kAna2020FHCPOT);
95 
96  // create experiment out of prediction and mock data
97  std::vector<const IChiSqExperiment*> expts;
98  expts.push_back(new SingleSampleExperiment(pred_numu,
99  mock_data));
100 
101  // include constrain on th13 from other experiments
102  expts.push_back(WorldReactorConstraint2019());
103 
104  // combine into one "MultiExperiment"
105  MultiExperiment * multi_experiment = new MultiExperiment(expts);
106 
107  // setup fitters
108  // fit vars for the constrained fit.
109  // does not include th23 because our null hypothesis
110  // specifies this value and so it is fixed.
111  std::vector<const IFitVar*> constrained_fit_vars;
112  constrained_fit_vars.push_back(&kFitDmSq32);
113  constrained_fit_vars.push_back(&kFitSinSq2Theta13);
114  MinuitFitter constrained_fit(multi_experiment,
115  constrained_fit_vars);
116 
117  // include th23 in the best fit to find the global maximum likelihood
118  std::vector<const IFitVar*> global_fit_vars;
119  global_fit_vars.push_back(&kFitSinSqTheta23);
120  global_fit_vars.push_back(&kFitDmSq32);
121  global_fit_vars.push_back(&kFitSinSq2Theta13);
122  // we may want to seed parameters at multiple points
123  // to avoid falling in local minima
124  SeedList global_fit_seeds({{&kFitSinSqTheta23, {0.4, 0.5, 0.6}}});
125  MinuitFitter global_fit(multi_experiment,
126  global_fit_vars);
127 
128 
130  SystShifts junk_shifts = SystShifts();
131  // run constrained fit
132  double chisq_constrained = constrained_fit.Fit(calc, IFitter::kQuiet);
133 
134  // save the results
135  double cf_dmsq = kFitDmSq32 .GetValue(calc);
136  double cf_delta = kFitDeltaInPiUnits.GetValue(calc);
137  double cf_th13 = kFitSinSq2Theta13 .GetValue(calc);
138 
140  // run global fit
141  double chisq_global = global_fit.Fit(calc,
142  junk_shifts,
143  global_fit_seeds,
144  {},
146  // save the results
147  double gf_dmsq = kFitDmSq32 .GetValue(calc);
148  double gf_delta = kFitDeltaInPiUnits.GetValue(calc);
149  double gf_th13 = kFitSinSq2Theta13 .GetValue(calc);
150  double gf_th23 = kFitSinSqTheta23 .GetValue(calc);
151 
152  FCPoint fcpt(delta, th23, th13, dmsq,
153  0, 0,
154  cf_delta, th23 , cf_th13, cf_dmsq, chisq_constrained,
155  gf_delta, gf_th23, gf_th13, gf_dmsq, chisq_global,
156  0, 0);
157 
158  fccol.AddPoint(fcpt);
159  }
160  std::cout << std::endl;
161  }
162 
163  fccol.SaveToFile(TString::Format("FCTutorial2020_FCCollection_%d_%d.root", jobid, procid).Data());
164 
165 }
166 
const XML_Char * name
Definition: expat.h:151
void PlotPseudoexperiments(std::vector< TH1 * > pse, std::string name)
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
Loads shifted spectra from files.
double GetValue(const osc::IOscCalcAdjustable *osc) const override
double th23
Definition: runWimpSim.h:98
double delta
Definition: runWimpSim.h:98
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Forward to wrapped Var&#39;s GetValue()
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:18
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Definition: FitVars.cxx:177
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Definition: FitVars.cxx:171
void FCTutorial2020(int NPSE=10, int nbins=1, int jobid=0, int procid=0)
void AddPoint(const FCPoint &p)
Definition: FCCollection.h:17
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Definition: FitVars.cxx:48
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Definition: FitVars.cxx:42
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
Spectrum MockData(double pot, int seed=0) const
Mock data is FakeData with Poisson fluctuations applied.
Definition: Spectrum.cxx:300
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
Represents the results of a single Feldman-Cousins pseudo-experiment.
Definition: FCPoint.h:8
const int nbins
Definition: cellShifts.C:15
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
const double kAna2020FHCPOT
Definition: Exposures.h:233
OStream cout
Definition: OStream.cxx:6
Combine multiple component experiments.
void SetProgress(double frac)
Update the progress fraction between zero and one.
Definition: Progress.cxx:41
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
A simple ascii-art progress bar.
Definition: Progress.h:9
const ReactorExperiment * WorldReactorConstraint2019()
Reactor constraint from PDG 2019.
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
void SaveToFile(const std::string &fname) const
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
double th13
Definition: runWimpSim.h:98
Compare a single data spectrum to the MC + cosmics expectation.
Collection of FCPoint. Serializable to/from a file.
Definition: FCCollection.h:12
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17
enum BeamMode string