run_joint_fit_2020_slices.C
Go to the documentation of this file.
3 #include "CAFAna/Fit/Fit.h"
11 #include "CAFAna/FC/FCSurface.h"
16 #include "CAFAna/Vars/FitVars.h"
17 #include "OscLib/IOscCalc.h"
18 #include "Utilities/rootlogon.C"
19 
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 void MakeMaps ( std::vector <const IFitVar*> profvars,
45  std::map<const IFitVar*, TGraph*> &profMap,
46  std::vector <const ISyst* > systs,
47  std::map<const ISyst*, TGraph*> &systMap);
48 void SaveMaps ( TDirectory * dir,
49  TString prefix,
50  std::vector <const IFitVar*> profvars,
51  std::vector <TString > profnames,
52  std::map<const IFitVar*, TGraph*> &profMap,
53  std::vector <const ISyst* > systs,
54  std::map<const ISyst*, TGraph*> &systMap);
55 
56 void run_joint_fit_2020_slices(bool corrSysts = false,
57  bool runOnGrid = false,
58  TString options="joint_fake2019_both",
59  bool th23slice = false,
60  bool octantSlice = true,
61  bool dmsqSlice = false,
62  bool th13Slice = false)
63 
64 {
65  bool nueOnly = options.Contains("nueOnly");
66  bool numuOnly = options.Contains("numuOnly");
67  bool joint = options.Contains("joint");
68  assert (nueOnly || numuOnly || joint);
69 
70  bool FHCOnly = options.Contains("FHCOnly");
71  bool RHCOnly = options.Contains("RHCOnly");
72  bool both = options.Contains("both");
73  assert (FHCOnly || RHCOnly || both);
74 
75  bool fake2019 = options.Contains("fake2019");
76  bool realData = options.Contains("realData");
77 
78  auto suffix = options;
79  if(corrSysts) suffix+="_systs";
80  else suffix+="_stats";
81 
82  TString tag;
83  if(FHCOnly) tag="FHCOnly/";
84  if(RHCOnly) tag="RHCOnly/";
85  if(both) tag = "RHC_and_FHC/";
86 
87  if(th23slice) suffix+="_th23";
88  if(dmsqSlice) suffix+="_dmsq";
89  if(th13Slice) suffix+="_th13";
90  if(!octantSlice)suffix+="_noOct";
91 
92  TString outdir = "/nova/ana/nu_e_ana/Ana2020/prefits/";//"/nova/ana/nu_e_ana/Ana2020/Results/"+tag+"contours/FCInputs/";
93  if (runOnGrid) outdir = "./";
94 
95  TString outFCfilename (outdir + "slices_FCInput_2020_" + suffix);
96  if(joint && corrSysts && !th23slice && ! dmsqSlice && !th13Slice) outFCfilename = outdir + "slices_FCInput_2020_" +suffix;
97 
98  // outdir = "/nova/ana/nu_e_ana/Ana2020/Results/"+tag+"slices/";
99  //outdir += corrSysts?"syst/":"stat/";
100  //if (runOnGrid) outdir = "./";
101 
102  TString outfilename (outdir + "hist_slices_2020_" + suffix);
103 
104 
105  //////////////////////////////////////////////////
106  // Load Nue and Numu experiments
107  //////////////////////////////////////////////////
108  //need numu only for prestage seeds
109 
110  //std::string decomp = "noPt"; // make Pt extrap hardcoded
111  std::string decomp = ""; // make Pt extrap hardcoded
112  std::vector <ana::predictions> preds = LoadPredictions (corrSysts, runOnGrid, decomp, nueOnly || joint, numuOnly || joint, FHCOnly || both, RHCOnly || both);
113  std::vector <ana::predictions> preds_numu_only = LoadPredictions (corrSysts, runOnGrid, decomp, false, numuOnly || joint, FHCOnly || both, RHCOnly || both);
114  std::vector <Spectrum * > data = LoadRealData (runOnGrid, nueOnly || joint, numuOnly || joint, FHCOnly || both, RHCOnly || both);
115  std::vector <Spectrum * > data_numu_only = LoadRealData (runOnGrid, false, numuOnly || joint, FHCOnly || both, RHCOnly || both);
116  std::vector <const IExperiment * > expts;
117  std::vector <const IExperiment * > expts_numu_only;
118 
119  int nnumu = 4;
120 
121  auto calc_fake = DefaultOscCalc();
122  if(fake2019) SetFakeCalc(calc_fake, 0.565, 2.48e-3, 0);
123  else if(!realData) {std::cerr << "need setting for data\n"; exit(1);}
124 
125  for(int i = 0; i < int(preds.size()); ++i){
126  double POT = preds[i].pot;
127  auto thisdata = GetFakeData (preds[i].pred, calc_fake, POT, preds[i].cos.first, preds[i].livetime);
128  if(realData) thisdata = data[i];
129  auto hcos = preds[i].cos.first->ToTH1(preds[i].livetime, kBlack, kSolid, kLivetime);
130  cout<<i<<" "<<preds[i].name<<" POT "<<POT<<" tot MC "<<preds[i].pred->Predict(calc_fake).Integral(POT)<<
131  " cos "<<hcos->Integral()<<" cos er "<<preds[i].cos.second<<" analyze data "<<thisdata->Integral(thisdata->POT())<<endl;
132  expts.push_back(new SingleSampleExperiment(preds[i].pred, *thisdata, *preds[i].cos.first, preds[i].cos.second, true));
133  }
134 
135  cout<<"Make numu only experiment to get the seeds"<<endl;
136  //make numu only experiment for seeds:
137 
138  for(int i = 0; i < (int) preds_numu_only.size(); ++i){
139  double POT = preds_numu_only[i].pot;
140  auto thisdata = GetFakeData (preds_numu_only[i].pred, calc_fake, POT, preds_numu_only[i].cos.first, preds_numu_only[i].livetime);
141  if(realData) thisdata = data_numu_only[i];
142  auto hcos = preds_numu_only[i].cos.first->ToTH1(preds_numu_only[i].livetime, kBlack, kSolid, kLivetime);
143  cout<<i<<" "<<preds_numu_only[i].name<<" POT "<<POT<<" tot MC "<<preds_numu_only[i].pred->Predict(calc_fake).Integral(POT)<<
144  " cos "<<hcos->Integral()<<" cos er "<<preds_numu_only[i].cos.second<<" analyze data "<<thisdata->Integral(thisdata->POT())<<endl;
145  expts_numu_only.push_back(new SingleSampleExperiment(preds_numu_only[i].pred, *thisdata, *preds_numu_only[i].cos.first, preds_numu_only[i].cos.second, true));
146  }
147 
148  ////////////////////////////////////////////////////////////
149  // Add constraints, make experiments
150  ////////////////////////////////////////////////////////////
151  std::cout << "\nCreating multiexperiment\n" << std::endl;
152 
153  if (!th13Slice) {expts.push_back(WorldReactorConstraint2019()); std::cout << "Adding WorldReactorConstraint2019\n";}
154  else {std::cout << "No reactor constraint, that's th13 slice\n";}
155  if(nueOnly) {
156  std::cout << "Adding Dmsq32ConstraintPDG2019\n";
157  expts.push_back(&kDmsq32ConstraintPDG2019);
158  }
159  std::cout << "Creating Multiexperiment with a total of "
160  << expts.size() << " experiments\n\n" << std::endl;
161  auto exptThis = new MultiExperiment(expts);
162 
163  std::cout << "Creating Multiexperiment of numu only SimpleExp with a total of "
164  << expts_numu_only.size() << " experiments\n\n" << std::endl;
165  auto exptThis_numu_only = new MultiExperiment(expts_numu_only);
166 
167  ////////////////////////////////////////////////////////////
168  // Systematics
169  ////////////////////////////////////////////////////////////
170 
171  std::cout << "Systematics for the fit:\n";
172  bool ptExtrap = true;
173  if (decomp == "noPt") ptExtrap = false;
174  auto systs = GetJointFitSystematicList(corrSysts, nueOnly, numuOnly, FHCOnly||both, RHCOnly||both, true, ptExtrap);
175 
176  if(!nueOnly && ! numuOnly && corrSysts){
177  if(FHCOnly){
178  exptThis->SetSystCorrelations(0, GetCorrelations(true, true, ptExtrap));
179  auto notfornumu = GetCorrelations(false, true, ptExtrap);
180  for(int i =0; i < nnumu; ++i) exptThis->SetSystCorrelations(i+1, notfornumu);
181  }
182  if(RHCOnly){
183  exptThis->SetSystCorrelations(0, GetCorrelations(true, false, ptExtrap));
184  auto notfornumu = GetCorrelations(false, false, ptExtrap);
185  for(int i =0; i < nnumu; ++i) exptThis->SetSystCorrelations(i+1, notfornumu);
186  }
187  if(both){
188  exptThis->SetSystCorrelations(0, GetCorrelations(true, true, ptExtrap));
189  exptThis->SetSystCorrelations(1, GetCorrelations(true, false, ptExtrap));
190  auto notfornumufhc = GetCorrelations(false, true, ptExtrap);
191  for(int i =0; i < 4; ++i) exptThis->SetSystCorrelations(i+2, notfornumufhc);
192  auto notfornumurhc = GetCorrelations(false, false, ptExtrap);
193  for(int i =0; i < 4; ++i) exptThis->SetSystCorrelations(i+6, notfornumurhc);
194  }
195  }
196  if (nueOnly && corrSysts){
197  if (FHCOnly) exptThis->SetSystCorrelations(0, GetCorrelations(true, true, ptExtrap));
198  if (RHCOnly) exptThis->SetSystCorrelations(0, GetCorrelations(true, false, ptExtrap));
199  if (both) {
200  exptThis->SetSystCorrelations(0, GetCorrelations(true, true, ptExtrap));
201  exptThis->SetSystCorrelations(1, GetCorrelations(true, false, ptExtrap));
202  }
203  }
204 
205  std::cout << "Systematics for the numu only fit:\n";
206  auto systs_numu_only = GetJointFitSystematicList(corrSysts, false, true, FHCOnly||both, RHCOnly||both, true, ptExtrap);
207 
208  if (corrSysts && numuOnly && both){
209  auto notfornumufhc = GetCorrelations(false, true, ptExtrap);
210  for(int i =0; i < 4; ++i) {
211  exptThis->SetSystCorrelations(i, notfornumufhc);
212  exptThis_numu_only->SetSystCorrelations(i, notfornumufhc);
213  }
214  auto notfornumurhc = GetCorrelations(false, false, ptExtrap);
215  for(int i =0; i < 4; ++i) {
216  exptThis->SetSystCorrelations(i+4, notfornumurhc);
217  exptThis_numu_only->SetSystCorrelations(i+4, notfornumurhc);
218  }
219  }
220  ////////////////////////////////////////////////////////////
221  // Fit
222  ////////////////////////////////////////////////////////////
223  std::cout << "Starting the fit" << std::endl;
224 
226 
227  SystShifts auxShifts = SystShifts::Nominal();
228 
229  std::vector <SystShifts> seedShifts = {};
230  if (corrSysts){
231  for (double systshift:{-2, +2}){
232  SystShifts tempShifts (&kAnaLightlevelNDSyst, systshift); // attempt to fix flip between positive and negative pulls
233  seedShifts.push_back(tempShifts);
234  }
235  }
236  //////////////////////////////////////////////////////////////////////
237  ///////////////////////// Seed controller ////////////////////////////
238  //////////////////////////////////////////////////////////////////////
239 
240  cout<<"------------------- Start prestage seeds --------------------------"<<endl;
241 
242  double minchi_numu = 1E20;
243  double pre_seed_th23;
244  double pre_seed_dmsq;
245  ResetOscCalcToDefault(calc);
246  auxShifts.ResetToNominal();
247 
248  double maxmix = 0.514; // from the numu results
249  double numu_pre_seedLONH, numu_pre_seedUONH, numu_pre_seedLOIH, numu_pre_seedUOIH, dmsq_numu_pre_seedNH, dmsq_numu_pre_seedIH;
250 
251  for(int hie:{1}){
252  std::cout << "\n\nFinding test best fit " << (hie>0? "NH " : "IH ") << "\n\n";
253  MinuitFitter fitnumu_only(exptThis_numu_only, {&kFitSinSqTheta23, &kFitDmSq32}, {});
254 
255  auto thisminchi = fitnumu_only.Fit(calc, auxShifts ,
256  SeedList({{&kFitDmSq32,{hie*2.5e-3}},
257  {&kFitSinSqTheta23, {0.4} }}), {}, //seedShifts,
258  IFitter::kVerbose)->EvalMetricVal();
259 
260  if(thisminchi < minchi_numu) minchi_numu = thisminchi;
261 
262  pre_seed_th23 = kFitSinSqTheta23.GetValue(calc);
263  pre_seed_dmsq = kFitDmSq32.GetValue(calc);
264  }
265 
266  numu_pre_seedLONH = ((pre_seed_th23>maxmix)?(2*maxmix-pre_seed_th23):pre_seed_th23);
267  numu_pre_seedUONH = ((pre_seed_th23>maxmix)?pre_seed_th23:(2*maxmix-pre_seed_th23));
268 
269  ResetOscCalcToDefault(calc);
270  auxShifts.ResetToNominal();
271 
272  cout<<"------------------- End prestage seeds --------------------------"<<endl;
273 
274 
275  struct th23helper{
276  std::vector<double> seeds;
277  const IFitVar * var;
278  TString label;
279  };
280 
281 
282  std::vector <th23helper> th23seeds;
283 
284  th23helper temp;
285  temp.seeds = {0.3};
286  temp.var = &kFitSinSqTheta23LowerOctant;
287  temp.label = "LO";
288 
289  //for NH:
290  if(octantSlice) {
291  th23seeds.push_back( { {0.499, numu_pre_seedLONH}, &kFitSinSqTheta23LowerOctant, "LO"});
292  th23seeds.push_back( { {0.501, numu_pre_seedUONH}, &kFitSinSqTheta23UpperOctant, "UO"});
293  }
294  else th23seeds.push_back({ {numu_pre_seedLONH, 0.5, numu_pre_seedUONH}, &kFitSinSqTheta23, ""});
295 
296  std::vector<double> delta_seeds = {0, 0.5, 1., 1.5};
297  std::vector<double> th23_old_seeds = {0.45, 0.5, 0.55};
298 
299  /////////////////////////////////////////////////////////////////////////
300  ////////////////////////// Different best fit searches //////////////////
301  /////////////////////////////////////////////////////////////////////////
302 
303  // Find the best fit points
304  double minchi23 = 1E20;
305  for(int hie:{-1, 1}){
306  for (auto & thseed:th23seeds){
307 
308  std::cout << "\n\nFinding best fit " << (hie > 0 ? "NH " : "IH ")
309  << thseed.label << ", "
310  << "ssth23 seeds ";
311  for (auto const & s:thseed.seeds) std::cout << s << ", ";
312  std::cout << std::endl;
313 
314  std::vector <const IFitVar*> fitvars = {&kFitDeltaInPiUnits,
315  thseed.var,
316  &kFitDmSq32,
318 
319  MinuitFitter fit23(exptThis, fitvars, systs);
320  cout<<" pre dmsq seed is "<<pre_seed_dmsq<<endl;
321  auto thisminchi = fit23.Fit(calc, auxShifts,
322  SeedList({
323  { &kFitDmSq32, {hie*pre_seed_dmsq} },
324  { thseed.var, thseed.seeds },
325  { &kFitDeltaInPiUnits, delta_seeds }
326  }), seedShifts,
327  IFitter::kVerbose)->EvalMetricVal();
328  minchi23= min(minchi23, thisminchi);
329 
330  if(corrSysts){
331  // Plot the systematic pulls and label with BFV
332  auto shifts = PlotSystShifts(auxShifts);
333  TString str = "Best fit " ;
334  for (auto &v:fitvars){
335  str += TString::Format(" %s=%.3f ",v->LatexName().c_str(),v->GetValue(calc));
336  }
337  str+= TString::Format(" LL=%.6f", thisminchi);
338  shifts->SetTitle(str);
339  gPad->Update();
340  TLine *l=new TLine(gPad->GetUxmin(),0,gPad->GetUxmax(),0);
341  l->Draw("same");
342  gPad->Print(outdir + "debug_slice_shifts_bestfit_" + suffix +
343  (hie>0? "NH": "IH") + thseed.label + ".pdf");
344  }
345 
346  ResetOscCalcToDefault(calc);
347  auxShifts.ResetToNominal();
348  }
349  }
350  std::cout << "\nFound overall minchi " << minchi23 << "\n\n";
351 
352  cout<<"Check with oldstyle seeds "<<endl;
353  double minchi23test = 1E20;
354  for(int hie:{-1,1}){
355  std::cout << "\n\nFinding test best fit "
356  << (hie>0? "NH " : "IH ")
357  << "\n\n";
358  std::vector <const IFitVar*> fitvars = {&kFitDeltaInPiUnits,
360  &kFitDmSq32,
362 
363  MinuitFitter fit23(exptThis, fitvars, systs);
364  auto thisminchi = fit23.Fit(calc,auxShifts ,SeedList({
365  {&kFitDmSq32,{hie*fabs(kFitDmSq32.GetValue(calc))}},
366  {&kFitSinSqTheta23, th23_old_seeds},
367  {&kFitDeltaInPiUnits, delta_seeds}
368  }), seedShifts
369  )->EvalMetricVal();
370  minchi23test= min(minchi23test, thisminchi);
371  if(corrSysts){
372  auto shifts = PlotSystShifts(auxShifts);
373  TString str = "Best fit " ;
374  for (auto &v:fitvars){
375  str += TString::Format(" %s=%.3f ",v->LatexName().c_str(),v->GetValue(calc));
376  }
377  str+= TString::Format(" LL=%.6f", thisminchi);
378  shifts->SetTitle(str);
379  gPad->Print(outdir + "debug_slice_shifts_bestfit_" + suffix +
380  (hie>0? "NH": "IH") + "test_chisq_with_no_octant_differentiation.pdf");
381  }
382 
383  ResetOscCalcToDefault(calc);
384  auxShifts.ResetToNominal();
385  }//end hie
386  std::cout << "\nFound overall test minchi " << minchi23test << std::endl;
387 
388  /////////////////////////////////////////////////////////////////////////////
389  /////////////////////// Now, slices! ////////////////////////////////////////
390  ////////////////////////////////////////////////////////////////////////////
391 
392  TFile * outfile = new TFile(outfilename+".root","recreate");
393  TFile * outFCfile = new TFile(outFCfilename+".root","recreate");
394 
395  outfile->cd();
396  TVectorD v(1);
397  v[0] = minchi23;
398  v.Write("overall_minchi");
399 
400  int steps = 60;// 80;
401 
402  //Now create all the slices
403  for(int hie: {-1, +1}){
404 
405  ResetOscCalcToDefault(calc);
406  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
407 
408  std::cout << "Starting profile " << (hie>0? "NH ": "IH") << "\n\n";
409 
410  struct ProfDef{
411  const IFitVar * fitvar;
412  double minX;
413  double maxX;
414  std::vector <const IFitVar * > profvars;
415  std::vector <TString> profvarnames;
416  std::map <const IFitVar *, std::vector <double> >profseeds;
417  TString shortname;
418  };
419 
420  for(auto const & thseed:th23seeds){
421  ProfDef profdef;
422  const IFitVar* kCurTheta23 = thseed.var;
423 
424  if(!th23slice && !dmsqSlice && !th13Slice){
425  profdef.fitvar = &kFitDeltaInPiUnits;
426  profdef.minX = 0;
427  profdef.maxX = 2;
428  profdef.profvars = {&kFitDmSq32, &kFitSinSq2Theta13, kCurTheta23};
429  profdef.profvarnames = {"DmSq32","SinSq2Theta13","SinSqTheta23"};
430  profdef.profseeds = {{kCurTheta23, thseed.seeds}, { &kFitDmSq32, {hie*pre_seed_dmsq}}, {&kFitSinSq2Theta13, {0.082}} };
431  profdef.shortname="delta";
432  }
433 
434 
435  if(th23slice){
436  profdef.fitvar = &kFitSinSqTheta23;
437  profdef.minX = (nueOnly ? 0:0.3);
438  profdef.maxX = (nueOnly ? 1:0.7);
439  profdef.profvars = {&kFitDmSq32, &kFitSinSq2Theta13, &kFitDeltaInPiUnits};
440  profdef.profvarnames = {"DmSq32","SinSq2Theta13","DeltaCPpi"};
441  profdef.profseeds = { {&kFitDeltaInPiUnits,delta_seeds}};
442  profdef.shortname="th23";
443  }
444 
445  if(th13Slice){
446  profdef.fitvar = &kFitSinSq2Theta13;
447  profdef.minX = 0.01;
448  profdef.maxX = 0.2;
449  profdef.profvars = {&kFitDmSq32, kCurTheta23, &kFitDeltaInPiUnits};
450  profdef.profvarnames = {"DmSq32","SinSqTheta23","DeltaCPpi"};
451  profdef.profseeds = { {kCurTheta23, thseed.seeds}, {&kFitDeltaInPiUnits,delta_seeds}, { &kFitDmSq32, {hie*pre_seed_dmsq}} };
452  profdef.shortname="th13";
453  }
454 
455  if(dmsqSlice){
456  profdef.fitvar = &kFitDmSq32;
457  profdef.minX = hie > 0 ? 2e-3 : -3e-3;
458  profdef.maxX = hie > 0 ? 3e-3 : -2e-3;
459  profdef.profvars = {kCurTheta23, &kFitSinSq2Theta13, &kFitDeltaInPiUnits};
460  profdef.profvarnames = {"SinSqTheta23","SinSq2Theta13","DeltaCPpi"};
461  profdef.profseeds = { {kCurTheta23, thseed.seeds}, {&kFitDeltaInPiUnits,delta_seeds}};
462  profdef.shortname="dmsq";
463  }
464 
465  std::map<const IFitVar*, TGraph*> profVarsMap;
466  std::map<const ISyst*, TGraph*> systMap;
467  MakeMaps (profdef.profvars, profVarsMap, systs, systMap);
468 
469  std::cout << " Profile "<< thseed.label << "\n";
470 
471  auto slice = SqrtProfile(exptThis, calc,
472  profdef.fitvar, steps, profdef.minX, profdef.maxX,
473  minchi23,
474  profdef.profvars,
475  systs,
476  profdef.profseeds,
477  seedShifts,
478  profVarsMap, systMap);
479  TString hieroctstr = (hie > 0 ? "NH" : "IH") + thseed.label;
480 
481  outfile->cd();
482  slice->Write((TString ) "slice_" + profdef.shortname + "_" + hieroctstr);
483  SaveMaps (outFCfile, hieroctstr,
484  profdef.profvars, profdef.profvarnames, profVarsMap, systs, systMap);
485 
486  }//end seeds
487  }//end hie
488 
489  delete outfile;
490  delete outFCfile;
491  std::cout << "\nProfiles saved to " << outfilename << std::endl;
492  return;
493 
494 }
495 
496 void MakeMaps ( std::vector <const IFitVar*> profvars,
497  std::map<const IFitVar*, TGraph*> &profMap,
498  std::vector <const ISyst* > systs,
499  std::map<const ISyst*, TGraph*> &systMap)
500 {
501 
502  for (const IFitVar * var:profvars)
503  profMap.insert(std::pair<const IFitVar*, TGraph*> (var, new TGraph()));
504  for (const ISyst* syst : systs)
505  systMap.insert(std::pair<const ISyst*, TGraph*> (syst, new TGraph()));
506 }
507 
508 void SaveMaps ( TDirectory * dir,
509  TString prefix,
510  std::vector <const IFitVar*> profvars,
511  std::vector <TString > profnames,
512  std::map<const IFitVar*, TGraph*> &profMap,
513  std::vector <const ISyst* > systs,
514  std::map<const ISyst*, TGraph*> &systMap)
515 {
516  TDirectory *tmp = gDirectory;
517  dir->cd();
518  for (int i = 0; i < (int) profvars.size(); ++i){
519  profMap[profvars[i]]->Write((prefix + "_" + profnames[i]));
520  }
521  for (const ISyst* s : systs)
522  systMap[s]->Write((prefix + "_" + s->ShortName()));
523  tmp->cd();
524 }
525 
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
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
double GetValue(const osc::IOscCalcAdjustable *osc) const override
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:18
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Definition: FitVars.cxx:171
const FitSinSqTheta23UpperOctant kFitSinSqTheta23UpperOctant
Definition: FitVars.cxx:16
void SaveMaps(TDirectory *out, std::map< std::string, IDecomp * > decomps_nominal, std::map< std::string, std::map< std::string, std::map< int, IDecomp * > > > decomps_shifted, std::map< std::string, PredictionNoExtrap * > predsNE_nominal, std::map< std::string, std::map< std::string, std::map< int, PredictionNoExtrap * > > > predsNE_shifted, std::map< std::string, PredictionSterile * > predsSt_nominal, std::map< std::string, std::map< std::string, std::map< int, PredictionSterile * > > > predsSt_shifted)
Save all of the objects in the input maps to the out directory/file.
Definition: PPFXHelper.h:1077
const Dmsq32Constraint kDmsq32ConstraintPDG2019(2.444e-3, 0.034e-3, 2.55e-3, 0.04e-3)
OStream cerr
Definition: OStream.cxx:7
std::vector< predictions > LoadPredictions(bool corrSysts=false, bool runOnGrid=false, std::string decomp="", bool nue=true, bool numu=true, bool fhc=true, bool rhc=true, bool NERSC=false)
Float_t tmp
Definition: plot.C:36
void MakeMaps(std::vector< const IFitVar * > profvars, std::map< const IFitVar *, TGraph * > &profMap, std::vector< const ISyst * > systs, std::map< const ISyst *, TGraph * > &systMap)
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
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
virtual void SetDmsq32(const T &dmsq32)=0
const char * label
const XML_Char const XML_Char * data
Definition: expat.h:268
const XML_Char * s
Definition: expat.h:262
string outfilename
knobs that need extra care
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:64
virtual T GetDmsq32() const
Definition: IOscCalc.h:86
void SetFakeCalc(osc::IOscCalcAdjustable *calc, double ssth23=-5, double dmsq32=-5, double dCP_Pi=-5)
const XML_Char * prefix
Definition: expat.h:380
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
std::vector< Spectrum * > LoadRealData(bool runOnGrid=false, bool nue=true, bool numu=true, bool fhc=true, bool rhc=true)
Combine multiple component experiments.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
TDirectory * dir
Definition: macro.C:5
double livetime
Definition: saveFDMCHists.C:21
const FitSinSqTheta23LowerOctant kFitSinSqTheta23LowerOctant
Definition: FitVars.cxx:17
void ResetToNominal()
Definition: SystShifts.cxx:141
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
exit(0)
T cos(T number)
Definition: d0nt_math.hpp:78
assert(nhit_max >=nhit_nbins)
Interface definition for fittable variables.
Definition: IFitVar.h:16
std::vector< std::pair< const ISyst *, const ISyst * > > GetCorrelations(bool isNue, bool isFHC, bool ptExtrap)
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
const DummyAnaSyst kAnaLightlevelNDSyst("Light_Level_ND","Light_Level_ND")
void run_joint_fit_2020_slices(bool corrSysts=false, bool runOnGrid=false, TString options="joint_fake2019_both", bool th23slice=false, bool octantSlice=true, bool dmsqSlice=false, bool th13Slice=false)
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
FILE * outfile
Definition: dump_event.C:13
Compare a single data spectrum to the MC + cosmics expectation.
TGraph * SqrtProfile(const IExperiment *expt, osc::IOscCalcAdjustable *calc, const IFitVar *v, int nbinsx, double minx, double maxx, double minchi, std::vector< const IFitVar * > profVars, std::vector< const ISyst * > profSysts, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts, std::map< const IFitVar *, TGraph * > &profVarsMap, std::map< const ISyst *, TGraph * > &systsMap, MinuitFitter::FitOpts opts)
Forward to Profile but sqrt the result for a crude significance.
Definition: Fit.cxx:143
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17