Tutorial2019Fit.C
Go to the documentation of this file.
3 #include "CAFAna/Analysis/Fit.h"
5 #include "CAFAna/Analysis/NuePlotStyle.h"
6 #include "CAFAna/Analysis/Surface.h"
12 #include "CAFAna/FC/FCSurface.h"
13 #include "CAFAna/Prediction/PredictionSystJoint2018.h"
15 #include "CAFAna/Systs/NumuSysts.h"
16 #include "CAFAna/Vars/FitVars.h"
17 #include "OscLib/IOscCalc.h"
18 #include "Utilities/rootlogon.C"
19 
20 #include "CAFAna/nue/Ana2019/joint_fit_2019_loader_tools.h"
21 
22 #include "TCanvas.h"
23 #include "TBox.h"
24 #include "TColor.h"
25 #include "TGraph.h"
26 #include "TVectorD.h"
27 #include "TF1.h"
28 #include "TLegend.h"
29 #include "TText.h"
30 #include "TLatex.h"
31 #include "TPad.h"
32 #include "TLine.h"
33 #include "TMarker.h"
34 #include "TStyle.h"
35 #include "TSystem.h"
36 #include "TGaxis.h"
37 
38 #include <algorithm>
39 #include <vector>
40 #include <string>
41 
42 using namespace ana;
43 
44 std::vector <Spectrum * > GetMockData(TString ana_options);
45 
46 void Tutorial2019Fit(bool corrSysts = true, /// do you want a fit with systematics or not?
47  TString options="joint_realData_both") /// what kind of fit is it:
48  /// joint | nueOnly | numuOnly
49  /// with loading realData or creating fake2019 data?
50  /// with both | RHCOnly | FHCOnly mode(s)?
51 
52 {
53  /// some logic that will be used further
54  bool nueOnly = options.Contains("nueOnly");
55  bool numuOnly = options.Contains("numuOnly");
56  bool joint = options.Contains("joint");
57  assert (nueOnly || numuOnly || joint);
58 
59  bool FHCOnly = options.Contains("FHCOnly");
60  bool RHCOnly = options.Contains("RHCOnly");
61  bool both = options.Contains("both");
62  assert (FHCOnly || RHCOnly || both);
63 
64  bool fake2019 = options.Contains("fake2019");
65  bool realData = options.Contains("realData");
66 
67  TString outdir = "./";
68 
69  //////////////////////////////////////////////////
70  // Load Nue and Numu experiments
71  //////////////////////////////////////////////////
72 
73  struct predictions {
74  const string name;
75  const IPrediction * pred;
76  std::pair <TH1D*, double> cos;
77  double pot;
78  };
79 
80  /// for making SingleExperiment we will need predictions (MC and cosmic) + data
81 
82  std::vector <predictions> preds;
83  std::vector <Spectrum * > data;
84  std::vector <const IExperiment * > expts;
85 
86  /// in case of making fake data set a calculator
87  auto calc_fake = DefaultOscCalc();
88  if(fake2019) SetFakeCalc(calc_fake, 0.565, 2.48e-3, 1.99);
89 
90  /// in case of joint ana the order in the prediciton vector is the next:
91  /// 0 - nue FHC, 1 - nue RHC, [2, 5] - 4Q for numu FHC and [6,9] - 4Q for numu RHC
92 
93  /// form prediction vectors depending on the type of fit you chose
94  if(!numuOnly) {
95  if(FHCOnly || both ) {
96  preds.push_back({"nue_fhc", GetNuePrediction2019("combo", DefaultOscCalc(), corrSysts, "fhc", false), GetNueCosmics2019("fhc"), GetPOT()}); /// prediction downloading functions are defined in the ..._tools.h header, GetPOT () returns the POT exposure value for rhc or fhc mode
97  }
98  if(RHCOnly || both ) {
99  preds.push_back({"nue_fhc", GetNuePrediction2019("prop", DefaultOscCalc(), corrSysts, "rhc", false), GetNueCosmics2019("rhc"), GetPOT(false)});
100  }
101  }
102 
103  int nnumu = 4;
104  if(!nueOnly) {
105  if(FHCOnly || both ) {
106  auto numu_preds = GetNumuPredictions2019(nnumu, corrSysts, "fhc");
107  auto numu_cosmics = GetNumuCosmics2019(nnumu, "fhc");
108  for (int i = 0; i< nnumu; i++ ){
109  preds.push_back({"numu_fhc_"+std::to_string(i+1), numu_preds[i], numu_cosmics[i], GetPOT()});
110  }
111  }
112  if(RHCOnly || both ) {
113  auto numu_preds = GetNumuPredictions2019(nnumu, corrSysts, "rhc");
114  auto numu_cosmics = GetNumuCosmics2019(nnumu, "rhc");
115  for (int i = 0; i< nnumu; i++ ){
116  preds.push_back({"numu_rhc_"+std::to_string(i+1), numu_preds[i], numu_cosmics[i], GetPOT(false)});
117  }
118  }
119  }
120 
121  /// Load data for all predictions
122 
123  if(realData){
124  if(!numuOnly){
125  if(FHCOnly || both ) {
126  auto temp = GetMockData("nue_fhc");
127  data.insert(data.end(),temp.begin(), temp.end());
128  }
129  if(RHCOnly || both ) {
130  auto temp = GetMockData("nue_rhc");
131  data.insert(data.end(),temp.begin(), temp.end());
132  }
133  }
134  if(!nueOnly){
135  if(FHCOnly || both ){
136  auto temp = GetMockData("numu_fhc");
137  data.insert(data.end(), temp.begin(), temp.end());}
138  if(RHCOnly || both ){
139  auto temp = GetMockData("numu_rhc");
140  data.insert(data.end(), temp.begin(), temp.end());}
141  }
142  }
143 
144 
145  /// Lets make collection of SingleExperiments:
146 
147  for(int i = 0; i < int(preds.size()); ++i){
148  double POT = preds[i].pot;
149  auto thisdata = GetFakeData (preds[i].pred, calc_fake, POT, preds[i].cos.first); /// create fake data if you chose this option
150  if(realData) thisdata = data[i]; /// if "realData" then replace by realData
151  cout<<i<<" "<<preds[i].name<<" POT "<<POT<<" tot MC "<<preds[i].pred->Predict(calc_fake).Integral(POT)<<" cos "<<preds[i].cos.first->Integral()<<" cos er "<<preds[i].cos.second<<" analyze data "<<thisdata->Integral(POT)<<endl;
152  expts.push_back(new SingleSampleExperiment(preds[i].pred, *thisdata, preds[i].cos.first, preds[i].cos.second));
153  /// debug plot to visualize what predictions and data you put for each experiment
154  new TCanvas();
155  DataMCComparison (*thisdata, preds[i].pred, DefaultOscCalc());
156  preds[i].cos.first->SetMarkerStyle(34);
157  preds[i].cos.first->SetMarkerColor(kCosmicBackgroundColor);
158  preds[i].cos.first->Draw("hist p same");
159  gPad->Print(outdir + "debug_predictions_" + std::to_string(i) + ".pdf");
160  }
161 
162  ////////////////////////////////////////////////////////////
163  // Add constraints, make experiments
164  ////////////////////////////////////////////////////////////
165  /// we use reactor experiment constraint for the data fit - it adds one more pull term with PDG values
166  std::cout << "\nCreating multiexperiment\n" << std::endl;
167  expts.push_back(WorldReactorConstraint2017()); std::cout << "Adding WorldReactorConstraint2017\n";
168 
169  if(nueOnly) {
170  std::cout << "Adding Dmsq32ConstraintPDG2017\n"; /// if we don't use our own numu data then put a constraint for nueOnly fit with PDG values
171  expts.push_back(&kDmsq32ConstraintPDG2017);
172  }
173 
174  std::cout << "Creating Multiexperiment with a total of "
175  << expts.size() << " experiments\n\n" << std::endl;
176  auto exptThis = new MultiExperiment(expts);
177 
178  ////////////////////////////////////////////////////////////
179  // Systematics
180  ////////////////////////////////////////////////////////////
181  /// Lets follow analysis here as well:
182  std::cout << "Systematics for the fit:\n";
183  auto systs = GetJointFitSystematicList(corrSysts, nueOnly, numuOnly, FHCOnly||both, RHCOnly||both, true);
184 
185  std::cout << "\n\nSystematic correlations...\n";
186  if(!nueOnly && ! numuOnly && corrSysts){
187  if(FHCOnly){
188  exptThis->SetSystCorrelations(0, GetCorrelations(true, true));
189  auto notfornumu = GetCorrelations(false, true);
190  for(int i =0; i < nnumu; ++i) exptThis->SetSystCorrelations(i+1, notfornumu);
191  }
192  if(RHCOnly){
193  exptThis->SetSystCorrelations(0, GetCorrelations(true, false));
194  auto notfornumu = GetCorrelations(false, true);
195  for(int i =0; i < nnumu; ++i) exptThis->SetSystCorrelations(i+1, notfornumu);
196  }
197  if(both){
198  exptThis->SetSystCorrelations(0, GetCorrelations(true, true));
199  exptThis->SetSystCorrelations(1, GetCorrelations(true, false));
200  auto notfornumu = GetCorrelations(false, true);
201  for(int i =0; i < 8; ++i) exptThis->SetSystCorrelations(i+2, notfornumu);
202  }
203  }
204  if (nueOnly && corrSysts){
205  if (FHCOnly) exptThis->SetSystCorrelations(0, GetCorrelations(true, true));
206  if (RHCOnly) exptThis->SetSystCorrelations(0, GetCorrelations(true, false));
207  if (both) {
208  exptThis->SetSystCorrelations(0, GetCorrelations(true, true));
209  exptThis->SetSystCorrelations(1, GetCorrelations(true, false));
210  }
211  }
212 
213  /// ... but it's too complicated and long for the tutorial, lets replace it by simplier and clear:
214  systs = {}; // empty vector in case stat. fit
215  if (corrSysts) systs = {&kAnaCalibrationSyst}; /// leave just one systematic here for syst. case
216 
217  ////////////////////////////////////////////////////////////
218  // Fit
219  ////////////////////////////////////////////////////////////
220  /// just preparation, reset all components
221  std::cout << "Starting the fit" << std::endl;
222 
224 
225  SystShifts auxShifts = SystShifts::Nominal();
226 
227  // In case some systs need different seeds
228  std::vector <SystShifts> seedShifts = {};
229 
230  //////////////////////////////////////////////////////////////////////
231  ///////////////////////// Seed controller ////////////////////////////
232  //////////////////////////////////////////////////////////////////////
233  /// one place to keep all seeds, a seed is a place where fitter will start work,
234  /// more seeds - more fit iterations and more chances to find a proper minima.
235  /// They help a lot to manage with difficult parameters. And a big drawback - they make fit longer.
236  /// So use it in smart way
237 
238  struct th23helper{
239  std::vector<double> seeds;
240  const IFitVar * var;
241  TString label;
242  };
243 
244  std::vector <th23helper> th23seeds;
245 
246  th23helper temp;
247  temp.seeds = {0.3};
248  temp.var = kFitSinSqTheta23LowerOctant;
249  temp.label = "LO";
250 
251  th23seeds.push_back( { {0.499, 0.45}, kFitSinSqTheta23LowerOctant, "LO"});
252  th23seeds.push_back( { {0.501, 0.55}, kFitSinSqTheta23UpperOctant, "UO"});
253 
254  double seed_dmsq = 2.5e-3;
255 
256  std::vector<double> delta_seeds = {0, 0.5, 1., 1.5};
257 
258  /////////////////////////////////////////////////////////////////////////
259  ////////////////////////// Different best fit searches //////////////////
260  /////////////////////////////////////////////////////////////////////////
261 
262  /// Find the best fit points
263  /// In NOvA we use the next scheme:
264  /// run {IH, NH} x {UO, LO} fits, find local minima for each one. Then choose the best fit (smallest LL)
265 
266  double minchi23 = 1E20;
267  for(int hie:{-1, 1}){
268  for (auto & thseed:th23seeds){
269 
270  std::cout << "\n\nFinding best fit " << (hie > 0 ? "NH " : "IH ")
271  << thseed.label << ", "
272  << "ssth23 seeds ";
273  for (auto const & s:thseed.seeds) std::cout << s << ", ";
274  std::cout << std::endl;
275  /// what parameters to vary:
276  std::vector <const IFitVar*> fitvars = {&kFitDeltaInPiUnits,
277  thseed.var, /// This is octant variable constraint in LO or UO
278  &kFitDmSq32,
280 
281  Fitter fit23(exptThis, fitvars, systs); /// define the fitter - what (multi-)experiment, what vars to vary, systematics
282  auto thisminchi = fit23.Fit(calc, auxShifts, /// do the job, define calulator (place where the BF values will be saved) and syst. shift vector
283  {
284  { &kFitDmSq32, {hie*seed_dmsq} }, /// seeds for specific variables
285  { thseed.var, thseed.seeds },
286  { &kFitDeltaInPiUnits, delta_seeds },
287  { &kFitSinSq2Theta13, {0.082}}
288  }, seedShifts, /// seeds for systematics
289  Fitter::kVerbose); /// kVerbose or kQuiet - detlization of minimization process
290  minchi23= min(minchi23, thisminchi); /// is it smaller or not than minchi23 which will contain the best chi
291 
292  if(corrSysts){
293  // Plot the systematic pulls and label with BFV
294  auto shifts = PlotSystShifts(auxShifts);
295  TString str = "Best fit " ;
296  for (auto &v:fitvars){
297  str += TString::Format(" %s=%.3f ",v->LatexName().c_str(),v->GetValue(calc));
298  }
299  str+= TString::Format(" LL=%.6f", thisminchi);
300  shifts->SetTitle(str);
301  gPad->Update();
302  TLine *l=new TLine(gPad->GetUxmin(),0,gPad->GetUxmax(),0);
303  l->Draw("same");
304  gPad->Print(outdir + "debug_slice_shifts_bestfit_" +
305  (hie>0? "NH": "IH") + thseed.label + ".pdf");
306  }
307  TString dirname = (hie > 0 ? "NH_" : "IH_") + thseed.label;
308  ResetOscCalcToDefault(calc); /// reset calc and syst. shift vector for the next round
309  auxShifts.ResetToNominal();
310  }
311  }
312  std::cout << "\nFound overall minchi " << minchi23 << "\n\n";
313 
314 } //end of main part
315 
316 std::vector <Spectrum * > GetMockData(TString ana_options){
317 
318  cout<<"\n\n ============= You're loading Mock data from the file ============="<<endl;
319  std::vector <Spectrum *> spec;
320  std::string Dir = "./";
321  bool RHCmode = ana_options.Contains("rhc");
322  bool FHCmode = ana_options.Contains("fhc");
323  bool NueData = ana_options.Contains("nue");
324  bool NumuData = ana_options.Contains("numu");
325  if (NueData){
326  if (FHCmode) spec.push_back(LoadFromFile<Spectrum>(Dir+"mock_data_for_yn_tutorial.root", "data_exp_0").release());
327  if (RHCmode) spec.push_back(LoadFromFile<Spectrum>(Dir+"mock_data_for_yn_tutorial.root", "data_exp_1").release());
328  }
329  if (NumuData){
330  if (FHCmode){
331  for (int i = 2; i < 6; ++i) {
332  auto temp = LoadFromFile<Spectrum>(Dir+"mock_data_for_yn_tutorial.root", "data_exp_"+std::to_string(i)).release();
333  spec.push_back(temp);
334  }
335  }
336  if (RHCmode){
337  for (int i = 6; i < 10; ++i) {
338  auto temp = LoadFromFile<Spectrum>(Dir+"mock_data_for_yn_tutorial.root", "data_exp_"+std::to_string(i)).release();
339  spec.push_back(temp);
340  }
341  }
342 
343  }
344  return spec;
345 
346 }
347 
348 
349 
Spectrum * GetFakeData(const IPrediction *pred, osc::IOscCalc *calc, const double pot, TH1D *cosmics=0)
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
TH1 * PlotSystShifts(const SystShifts &shifts, bool sortName)
Definition: Plots.cxx:1495
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
const Color_t kCosmicBackgroundColor
Definition: Style.h:41
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:18
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
std::vector< const IPrediction * > GetNumuPredictions2019(const int nq=4, bool useSysts=true, std::string beam="fhc", bool GetFromUPS=false, ENu2018ExtrapType numuExtrap=kNuMu, bool minimizeMemory=false, bool NERSC=false)
const Dmsq32Constraint kDmsq32ConstraintPDG2017(2.45e-3, 0.05e-3, 2.52e-3, 0.05e-3)
const FitSinSqTheta23UpperOctant kFitSinSqTheta23UpperOctant
Definition: FitVars.cxx:16
void Tutorial2019Fit(bool corrSysts=true, TString options="joint_realData_both")
static SystShifts Nominal()
Definition: SystShifts.h:34
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
const IPrediction * GetNuePrediction2019(std::string decomp, osc::IOscCalc *calc, bool corrSysts, std::string beam, bool isFakeData, bool GetFromUPS=false, bool minimizeMemory=false, bool NERSC=false)
const char * label
const XML_Char const XML_Char * data
Definition: expat.h:268
const XML_Char * s
Definition: expat.h:262
std::vector< Spectrum * > GetMockData(TString ana_options)
const DummyAnaSyst kAnaCalibrationSyst("Calibration","AbsCalib")
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
void SetFakeCalc(osc::IOscCalcAdjustable *calc, double ssth23=-5, double dmsq32=-5, double dCP_Pi=-5)
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
Combine multiple component experiments.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::vector< std::pair< TH1D *, double > > GetNumuCosmics2019(const int nq=4, std::string beam="fhc", bool GetFromUPS=false, bool NERSC=false)
const FitSinSqTheta23LowerOctant kFitSinSqTheta23LowerOctant
Definition: FitVars.cxx:17
void ResetToNominal()
Definition: SystShifts.cxx:141
assert(nhit_max >=nhit_nbins)
Interface definition for fittable variables.
Definition: IFitVar.h:16
double GetPOT(bool isFHC=true)
std::vector< std::pair< const ISyst *, const ISyst * > > GetCorrelations(bool isNue, bool isFHC, bool ptExtrap)
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
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
TH1 * DataMCComparison(const Spectrum &data, const Spectrum &mc, EBinType bintype)
Definition: Plots.cxx:35
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
std::pair< Spectrum *, double > cos
Float_t e
Definition: plot.C:35
#define for
Definition: msvc_pragmas.h:3
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
std::vector< const ISyst * > GetJointFitSystematicList(bool corrSysts, bool nueExclusive=false, bool numuExclusive=false, bool isFHC=true, bool isRHC=true, bool intersection=true, bool ptExtrap=true)
const std::string outdir
std::pair< TH1D *, double > GetNueCosmics2019(std::string beam, bool GetFromUPS=false, bool NERSC=false)
Compare a single data spectrum to the MC + cosmics expectation.