joint_fit_future_bestfit_univ.C
Go to the documentation of this file.
3 #include "CAFAna/Fit/Fit.h"
11 #include "CAFAna/FC/FCSurface.h"
15 #include "CAFAna/Vars/FitVars.h"
16 #include "Utilities/rootlogon.C"
17 #include "OscLib/IOscCalc.h"
18 
20 
21 #include "TCanvas.h"
22 #include "TMarker.h"
23 #include "TBox.h"
24 #include "TLatex.h"
25 #include "TColor.h"
26 #include "TGraph.h"
27 #include "TVectorD.h"
28 #include "TF1.h"
29 #include "TLegend.h"
30 #include "TLine.h"
31 #include "TMarker.h"
32 #include "TStyle.h"
33 #include "TSystem.h"
34 #include "TGaxis.h"
35 
36 #include <algorithm>
37 #include <vector>
38 #include <string>
39 #include <fstream>
40 #include <sstream>
41 
42 using namespace ana;
43 // script for making 1d sensitivities e.g. MH or CPV
44 void joint_fit_future_bestfit_univ(int universes=2,
45  TString uni_tag = "",
46  double deltaCP = 0.18,
47  double dmsq32 = 2.51e-3,
48  double th23 = 0.58,
49  string sign = "plus",
50  double FHCexp = kAna2020FHCPOT,
51  double RHCexp = kAna2020RHCPOT,
52  double FHClive = kAna2020FHCLivetime,
53  double RHClive = kAna2020RHCLivetime,
54  bool corrSysts = false,
55  TString options="joint_realData_both_onlyIH",
56  TString plot_type = "delta",
57  bool asimov = false)
58 {
59 
60  cout<<"\nLOAD PREDICTIONS\n\n"<<endl;
61 
62  bool nueOnly = options.Contains("nueOnly");
63  bool numuOnly = options.Contains("numuOnly");
64  bool joint = options.Contains("joint");
65  assert (nueOnly || numuOnly || joint);
66 
67  bool FHCOnly = options.Contains("FHCOnly");
68  bool RHCOnly = options.Contains("RHCOnly");
69  bool both = options.Contains("both");
70  assert (FHCOnly || RHCOnly || both);
71 
72  bool hieIH = options.Contains("IH");
73  bool hieNH = options.Contains("NH");
74  bool isDelta = plot_type.Contains("delta");
75 
76  assert (hieIH || hieNH);
77 
78  string syst = (corrSysts?"syst_":"stat_");
79 
80  TString outfilename (syst + "bestfit_for_"+ options + "_deltacp_" + std::to_string(deltaCP) + "_" + std::to_string(universes) + "_" + uni_tag);
81 
82  TFile * outfile = new TFile(outfilename+".root","recreate");
83 
84  //////////////////////////////////////////////////
85  // Load Nue and Numu experiments
86  //////////////////////////////////////////////////
87 
88  std::string decomp = ""; // use Pt extrap all the times
89  std::vector <ana::predictions> preds = LoadPredictions (corrSysts, true, decomp, nueOnly || joint, numuOnly || joint, FHCOnly || both, RHCOnly || both);
90  std::vector <ana::predictions> preds_numu_only = LoadPredictions (corrSysts, true, decomp, false, numuOnly || joint, FHCOnly || both, RHCOnly || both);
91  std::vector <Spectrum > data;
92  std::vector <Spectrum > data_numu_only;
93  std::vector <const IExperiment * > expts;
94  std::vector <const IExperiment * > expts_numu_only;
95 
96  auto calc_fake = DefaultOscCalc();
97  if (sign =="minus") dmsq32 = -dmsq32;
98  SetFakeCalc(calc_fake, th23, dmsq32, deltaCP);
99 
100  outfile->cd();
101  TVectorD osc(3);
102  osc[0] = th23;
103  osc[1] = dmsq32;
104  osc[2] = deltaCP;
105  osc.Write("osc_params");
106 
107  for(int j = 0; j< universes; j++){
108  cout<<"\nclear expts and data vectors\n"<<endl;
109  expts.clear();
110  expts_numu_only.clear();
111  data.clear();
112  data_numu_only.clear();
113  // isn't set for individual fits (only joint nue+numu rhc+fhc)
114  if(!numuOnly){
115  if(FHCOnly || both ) data.push_back(GenerateFutureData(calc_fake, preds[0].pred, preds[0].cos.first, FHCexp, FHClive));
116  if(RHCOnly || both ) data.push_back(GenerateFutureData(calc_fake, preds[1].pred, preds[1].cos.first, RHCexp, RHClive));
117  }
118  if(!nueOnly){
119  if(FHCOnly || both ){
120  std::vector <Spectrum> numu_data;
121  numu_data.push_back(GenerateFutureData(calc_fake, preds[2].pred, preds[2].cos.first, FHCexp, FHClive));
122  numu_data.push_back(GenerateFutureData(calc_fake, preds[3].pred, preds[3].cos.first, FHCexp, FHClive));
123  numu_data.push_back(GenerateFutureData(calc_fake, preds[4].pred, preds[4].cos.first, FHCexp, FHClive));
124  numu_data.push_back(GenerateFutureData(calc_fake, preds[5].pred, preds[5].cos.first, FHCexp, FHClive));
125  data.insert(data.end(),numu_data.begin(), numu_data.end());
126  data_numu_only.insert(data_numu_only.end(),numu_data.begin(), numu_data.end());
127  }
128  if(RHCOnly || both ){
129  std::vector <Spectrum> numu_data;
130  numu_data.push_back(GenerateFutureData(calc_fake, preds[6].pred, preds[6].cos.first, RHCexp, RHClive));
131  numu_data.push_back(GenerateFutureData(calc_fake, preds[7].pred, preds[7].cos.first, RHCexp, RHClive));
132  numu_data.push_back(GenerateFutureData(calc_fake, preds[8].pred, preds[8].cos.first, RHCexp, RHClive));
133  numu_data.push_back(GenerateFutureData(calc_fake, preds[9].pred, preds[9].cos.first, RHCexp, RHClive));
134  data.insert(data.end(),numu_data.begin(), numu_data.end());
135  data_numu_only.insert(data_numu_only.end(),numu_data.begin(), numu_data.end());
136  }
137  }
138 
139  cout<<"check all vector sizes: "<<preds.size()<<" data "<<data.size()<<endl;
140 
141  TDirectory* d = outfile->mkdir(("universe_"+std::to_string(j)).c_str());
142  outfile->cd(("universe_"+std::to_string(j)).c_str());
143 
144  for(int i = 0; i < int(preds.size()); ++i){
145  auto thisdata = data[i];
146  double POT = thisdata.POT();
147  double LT = thisdata.Livetime();
148  thisdata.SaveTo(d, ("data_exp_"+std::to_string(i)).c_str());
149  auto hcos = preds[i].cos.first->ToTH1(LT, kBlack, kSolid, kLivetime);
150  cout<<i<<" "<<preds[i].name<<" POT "<<POT<<" tot MC "<<preds[i].pred->Predict(calc_fake).Integral(POT)<<" cos "<<hcos->Integral()<<" cos er "<<preds[i].cos.second<<" analyze data "<<thisdata.Integral(POT)<<endl;
151  expts.push_back(new SingleSampleExperiment(preds[i].pred, thisdata, *preds[i].cos.first, preds[i].cos.second));
152  }
153 
154  cout<<"Make numu only experiment to get the seeds"<<endl;
155  //make numu only experiment for seeds:
156 
157  for(int i = 0; i < (int) preds_numu_only.size(); ++i){
158  auto thisdata = data_numu_only[i];
159  double POT = thisdata.POT();
160  double LT = thisdata.Livetime();
161  auto hcos = preds_numu_only[i].cos.first->ToTH1(LT, kBlack, kSolid, kLivetime);
162  cout<<i<<" "<<preds_numu_only[i].name<<" POT "<<POT<<" tot MC "<<preds_numu_only[i].pred->Predict(calc_fake).Integral(POT)<<" cos "<<hcos->Integral()<<" cos er "<<preds_numu_only[i].cos.second<<" analyze data "<<thisdata.Integral(POT)<<endl;
163  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));
164  }
165 
166  //////////////////////////////////////////////////
167  // Add constraints, make experiments
168  //////////////////////////////////////////////////
169  std::cout << "\nCreating multiexperiment\n" << std::endl;
170  std::cout << "Adding WorldReactorConstraint2019\n";
171  expts.push_back(WorldReactorConstraint2019());
172 
173  if(nueOnly) {
174  std::cout << "Adding Dmsq32ConstraintPDG2017\n";
175  expts.push_back(&kDmsq32ConstraintPDG2017);
176  }
177  std::cout << "Creating Multiexperiment with a total of "
178  << expts.size() << " experiments\n\n" << std::endl;
179  auto exptThis = new MultiExperiment(expts);
180 
181  std::cout << "Creating Multiexperiment of numu only SimpleExp with a total of "
182  << expts_numu_only.size() << " experiments\n\n" << std::endl;
183  auto exptThis_numu_only = new MultiExperiment(expts_numu_only);
184 
185 
186  ////////////////////////////////////////////////////////////
187  // Systematics
188  ////////////////////////////////////////////////////////////
189 
190  bool ptExtrap = true;
191  auto systs = GetJointFitSystematicList(corrSysts, nueOnly, numuOnly, FHCOnly||both, RHCOnly||both, true, ptExtrap);
192 
193  int nnumu = 4;
194  std::cout << "\n\nSystematic correlations...\n";
195  if(!nueOnly && ! numuOnly && corrSysts){
196  if(FHCOnly){
197  exptThis->SetSystCorrelations(0, GetCorrelations(true, true, ptExtrap));
198  auto notfornumu = GetCorrelations(false, true, ptExtrap);
199  for(int i =0; i < nnumu; ++i) exptThis->SetSystCorrelations(i+1, notfornumu);
200  }
201  if(RHCOnly){
202  exptThis->SetSystCorrelations(0, GetCorrelations(true, false, ptExtrap));
203  auto notfornumu = GetCorrelations(false, false, ptExtrap);
204  for(int i =0; i < nnumu; ++i) exptThis->SetSystCorrelations(i+1, notfornumu);
205  }
206  if(both){
207  exptThis->SetSystCorrelations(0, GetCorrelations(true, true, ptExtrap));
208  exptThis->SetSystCorrelations(1, GetCorrelations(true, false, ptExtrap));
209  auto notfornumufhc = GetCorrelations(false, true, ptExtrap);
210  for(int i =0; i < 4; ++i) exptThis->SetSystCorrelations(i+2, notfornumufhc);
211  auto notfornumurhc = GetCorrelations(false, false, ptExtrap);
212  for(int i =0; i < 4; ++i) exptThis->SetSystCorrelations(i+6, notfornumurhc);
213  }
214  }
215  if (nueOnly && corrSysts){
216  if (FHCOnly) exptThis->SetSystCorrelations(0, GetCorrelations(true, true, ptExtrap));
217  if (RHCOnly) exptThis->SetSystCorrelations(0, GetCorrelations(true, false, ptExtrap));
218  if (both) {
219  exptThis->SetSystCorrelations(0, GetCorrelations(true, true, ptExtrap));
220  exptThis->SetSystCorrelations(1, GetCorrelations(true, false, ptExtrap));
221  }
222  }
223 
224  std::cout << "Systematics for the numu only fit:\n";
225  auto systs_numu_only = GetJointFitSystematicList(corrSysts, false, true, FHCOnly||both, RHCOnly||both, true, ptExtrap);
226 
227  if (corrSysts && numuOnly && both){
228  auto notfornumufhc = GetCorrelations(false, true, ptExtrap);
229  for(int i =0; i < 4; ++i) {
230  exptThis->SetSystCorrelations(i, notfornumufhc);
231  exptThis_numu_only->SetSystCorrelations(i, notfornumufhc);
232  }
233  auto notfornumurhc = GetCorrelations(false, false, ptExtrap);
234  for(int i =0; i < 4; ++i) {
235  exptThis->SetSystCorrelations(i+4, notfornumurhc);
236  exptThis_numu_only->SetSystCorrelations(i+4, notfornumurhc);
237  }
238  }
239 
240  ////////////////////////////////////////////////////////////
241  // Fit
242  ////////////////////////////////////////////////////////////
243  std::cout << "Starting the fit" << std::endl;
244 
246 
247  SystShifts auxShifts = SystShifts::Nominal();
248  // In case some systs need different seeds
249  std::vector <SystShifts> seedShifts = {};
250 
251  //////////////////////////////////////////////////////////////////////
252  /////////////////////////// Seed controller ////////////////////////////
253  ////////////////////////////////////////////////////////////////////////
254 
255 
256  cout<<"------------------- Start prestage seeds --------------------------"<<endl;
257 
258  double minchi_numu = 1E20;
259  double pre_seed_th23;
260  double pre_seed_dmsq;
261  ResetOscCalcToDefault(calc);
262  auxShifts.ResetToNominal();
263 
264  double maxmix = 0.514; // from the numu results
265  double numu_pre_seedLONH, numu_pre_seedUONH, numu_pre_seedLOIH, numu_pre_seedUOIH, dmsq_numu_pre_seedNH, dmsq_numu_pre_seedIH;
266 
267  for(int hie:{1}){
268  std::cout << "\n\nFinding test best fit " << (hie>0? "NH " : "IH ") << "\n\n";
269  MinuitFitter fitnumu_only(exptThis_numu_only, {&kFitSinSqTheta23, &kFitDmSq32Scaled}, {});
270 
271  auto thisminchi = fitnumu_only.Fit(calc, auxShifts ,
272  SeedList({{&kFitDmSq32Scaled,{hie*2.5}},
273  {&kFitSinSqTheta23, {0.4} }}), {},// seedShifts
275 
276  if(thisminchi < minchi_numu) minchi_numu = thisminchi;
277  }
278  pre_seed_th23 = kFitSinSqTheta23.GetValue(calc);
279  pre_seed_dmsq = kFitDmSq32Scaled.GetValue(calc);
280 
281  numu_pre_seedLONH = ((pre_seed_th23>maxmix)?(2*maxmix-pre_seed_th23):pre_seed_th23);
282  numu_pre_seedUONH = ((pre_seed_th23>maxmix)?pre_seed_th23:(2*maxmix-pre_seed_th23));
283 
284  ResetOscCalcToDefault(calc);
285  auxShifts.ResetToNominal();
286 
287  cout<<"------------------- End prestage seeds --------------------------"<<endl;
288 
289  struct th23helper{
290  std::vector<double> seeds;
291  const IFitVar * var;
292  TString label;
293  };
294 
295  std::vector <th23helper> th23seeds;
296 
297  th23helper temp;
298  temp.seeds = {0.3};
299  temp.var = &kFitSinSqTheta23LowerOctant;
300  temp.label = "LO";
301 
302  //for NH:
303  th23seeds.push_back( { {0.499, numu_pre_seedLONH}, &kFitSinSqTheta23LowerOctant, "LO"});
304  th23seeds.push_back( { {0.501, numu_pre_seedUONH}, &kFitSinSqTheta23UpperOctant, "UO"});
305  std::vector<double> th23_all_seeds_NH = {numu_pre_seedLONH, 0.5, numu_pre_seedUONH};
306 
307  std::vector<double> delta_seeds = {0, 0.5, 1., 1.5};
308  std::vector<double> th23_old_seeds = {0.45, 0.5, 0.55};
309 
310  ///////////////////////////////////////////////////////////////////////////
311  //////////////////////////// Best fit searches ////////////////////////////
312  ///////////////////////////////////////////////////////////////////////////
313 
314  double minchi23 = 1E20; double thisdmsq, thisth23, thisdcp, thisth13;
315  double curdmsq, curth23, curdcp, curth13;
316 
317  for(int hie:{-1,1}){
318  for (auto & thseed: th23seeds){
319  std::cout << "\n\nFinding best fit "
320  << thseed.label
321  << (hie>0? " NH " : " IH ")
322  << "\n\n";
323  for (auto const & s:thseed.seeds) std::cout << s << ", ";
324  std::cout << std::endl;
325  cout<<"pre seed from numu is "<<pre_seed_dmsq<<endl;
326  std::vector <const IFitVar*> fitvars = {&kFitDeltaInPiUnits,
327  thseed.var,
330 
331 
332  MinuitFitter fit23(exptThis, fitvars, systs);
333  auto thisminchi = fit23.Fit(calc,auxShifts ,SeedList({
334  {&kFitDmSq32Scaled,{hie*pre_seed_dmsq}},
335  {thseed.var, thseed.seeds},
336  {&kFitDeltaInPiUnits, delta_seeds},
337  {&kFitSinSq2Theta13, {calc->GetTh13()}} }),
339  );
340 
341  if(corrSysts){
342  auto shifts = PlotSystShifts(auxShifts);
343  TString str = "Best fit " ;
344  for (auto &v:fitvars){
345  str += TString::Format(" %s=%.3f ",v->LatexName().c_str(),v->GetValue(calc));
346  }
347  str+= TString::Format(" LL=%.6f", thisminchi);
348  shifts->SetTitle(str);
349  string hietag = (hie>0? "NH": "IH");
350  gPad->Print("debug_slice_shifts_bestfit_" + options + "_univ_" + std::to_string(j) + hietag + thseed.label + ".pdf");
351  }
352  curdmsq = kFitDmSq32Scaled.GetValue(calc);
353  curth13 = kFitSinSq2Theta13.GetValue(calc);
354  curth23 = thseed.var->GetValue(calc);
355  curdcp = kFitDeltaInPiUnits.GetValue(calc);
356  TVectorD v(5);
357  v[0] = thisminchi;
358  v[1] = curdmsq;
359  v[2] = curth23;
360  v[3] = curdcp;
361  v[4] = curth13;
362  string hietag = (hie>0? "NH": "IH");
363  v.Write("chi_"+std::to_string(j)+"_"+hietag+"_"+thseed.label);
364  if(thisminchi<minchi23){
365  thisdmsq = kFitDmSq32Scaled.GetValue(calc);
366  thisth13 = kFitSinSq2Theta13.GetValue(calc);
367  thisth23 = thseed.var->GetValue(calc);
368  thisdcp = kFitDeltaInPiUnits.GetValue(calc);
369  cout <<"this values (dmsq, th13, th23, dcp): "<<thisdmsq<<" "<<thisth13<<" "<<thisth23<<" "<<thisdcp<<endl;
370  }
371  minchi23= min(minchi23, thisminchi);
372  ResetOscCalcToDefault(calc);
373  auxShifts.ResetToNominal();
374  }//end th23
375  }//end hie
376  std::cout << "\nFound overall minchi " << minchi23 << std::endl;
377 
378  TVectorD v(5);
379  v[0] = minchi23;
380  v[1] = thisdmsq;
381  v[2] = thisth23;
382  v[3] = thisdcp;
383  v[4] = thisth13;
384  cout<<"minchi wrote into the file "<<minchi23<<endl;
385  v.Write(("bf_"+std::to_string(j)).c_str());
386 
387  if(isDelta){
388  double minchi23_delta = 1E20;
390  for(int hie:{-1,1}){
391  calc = DefaultOscCalc();
392  calc->SetdCP(deltaCP*M_PI);
393  auto minchi = fit23_dcp.Fit(calc, auxShifts , SeedList({
394  {&kFitDmSq32Scaled,{hie*pre_seed_dmsq}},
395  {&kFitDmSq32Scaled, th23_old_seeds},
396  {&kFitSinSq2Theta13, {calc->GetTh13()}} }),
398  );
399  minchi23_delta = min(minchi23_delta, minchi);
400  curdmsq = kFitDmSq32Scaled.GetValue(calc);
401  curth13 = kFitSinSq2Theta13.GetValue(calc);
402  curth23 = kFitSinSqTheta23.GetValue(calc);
403  curdcp = kFitDeltaInPiUnits.GetValue(calc);
404  TVectorD v(5);
405  v[0] = minchi;
406  v[1] = curdmsq;
407  v[2] = curth23;
408  v[3] = curdcp;
409  v[4] = curth13;
410  string hietag = (hie>0? "NH": "IH");
411  v.Write(("this_delta_chi_"+std::to_string(j)+"_"+hietag).c_str());
412  if(minchi<minchi23_delta){
413  thisdmsq = kFitDmSq32Scaled.GetValue(calc);
414  thisth13 = kFitSinSq2Theta13.GetValue(calc);
415  thisth23 = kFitSinSqTheta23.GetValue(calc);
416  thisdcp = kFitDeltaInPiUnits.GetValue(calc);
417  cout <<"this values (dmsq, th13, th23, dcp): "<<thisdmsq<<" "<<thisth13<<" "<<thisth23<<" "<<thisdcp<<endl;
418  }
419  ResetOscCalcToDefault(calc);
420  auxShifts.ResetToNominal();
421 
422  }
423  double minchi23_cpv = 1E20;
424  for(int hie:{-1,1}){
425  for(int dcp:{0,1}) {
426  calc = DefaultOscCalc();
427  calc->SetdCP(dcp*M_PI);
428  auto minchi = fit23_dcp.Fit(calc,auxShifts , SeedList({
429  {&kFitDmSq32Scaled,{hie*pre_seed_dmsq}},
430  {&kFitDmSq32Scaled, th23_old_seeds},
431  {&kFitSinSq2Theta13, {calc->GetTh13()}} }),
433  );
434  minchi23_cpv = min(minchi23_cpv, minchi);
435  curdmsq = kFitDmSq32Scaled.GetValue(calc);
436  curth13 = kFitSinSq2Theta13.GetValue(calc);
437  curth23 = kFitSinSqTheta23.GetValue(calc);
438  curdcp = kFitDeltaInPiUnits.GetValue(calc);
439  TVectorD v(5);
440  v[0] = minchi;
441  v[1] = curdmsq;
442  v[2] = curth23;
443  v[3] = curdcp;
444  v[4] = curth13;
445  string hietag = (hie>0? "NH": "IH");
446  string dcp_tag = (dcp>0? "pi": "zero");
447  v.Write(("cpv_chi_"+std::to_string(j)+"_"+hietag+"_"+dcp_tag).c_str());
448  if(minchi<minchi23_cpv){
449  thisdmsq = kFitDmSq32Scaled.GetValue(calc);
450  thisth13 = kFitSinSq2Theta13.GetValue(calc);
451  thisth23 = kFitSinSqTheta23.GetValue(calc);
452  thisdcp = kFitDeltaInPiUnits.GetValue(calc);
453  cout <<"this values (dmsq, th13, th23, dcp): "<<thisdmsq<<" "<<thisth13<<" "<<thisth23<<" "<<thisdcp<<endl;
454  }
455  ResetOscCalcToDefault(calc);
456  auxShifts.ResetToNominal();
457  }
458  }
459 
460  TVectorD v_cpv(5);
461  v_cpv[0] = minchi23_cpv;
462  v_cpv[1] = thisdmsq;
463  v_cpv[2] = thisth23;
464  v_cpv[3] = thisdcp;
465  v_cpv[4] = thisth13;
466  cout<<"cpv minchi wrote into the file "<<minchi23_cpv<<endl;
467  v_cpv.Write(("bf_cpv_"+std::to_string(j)).c_str());
468  }
469 
470  cout<<"\n\n done with universe no "<<j<<"\n\n\n";
471  }
472 }
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
double GetValue(const osc::IOscCalcAdjustable *osc) const override
double th23
Definition: runWimpSim.h:98
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Forward to wrapped Var&#39;s GetValue()
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
const Dmsq32Constraint kDmsq32ConstraintPDG2017(2.45e-3, 0.05e-3, 2.52e-3, 0.05e-3)
const FitSinSqTheta23UpperOctant kFitSinSqTheta23UpperOctant
Definition: FitVars.cxx:16
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)
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
Spectrum GenerateFutureData(osc::IOscCalc *calc, const IPrediction *pred, Spectrum *cosmics, double futurePOT, double futureLT, bool randomSystShifts=false, bool asimov=false)
#define M_PI
Definition: SbMath.h:34
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
const double kAna2020FHCLivetime
Definition: Exposures.h:237
Float_t d
Definition: plot.C:236
const double j
Definition: BetheBloch.cxx:29
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
const double kAna2020FHCPOT
Definition: Exposures.h:233
void SetFakeCalc(osc::IOscCalcAdjustable *calc, double ssth23=-5, double dmsq32=-5, double dCP_Pi=-5)
const double kAna2020RHCPOT
Definition: Exposures.h:235
Oscillation probability calculators.
Definition: Calcs.h: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
virtual T GetTh13() const
Definition: IOscCalc.h:88
Combine multiple component experiments.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void joint_fit_future_bestfit_univ(int universes=2, TString uni_tag="", double deltaCP=0.18, double dmsq32=2.51e-3, double th23=0.58, string sign="plus", double FHCexp=kAna2020FHCPOT, double RHCexp=kAna2020RHCPOT, double FHClive=kAna2020FHCLivetime, double RHClive=kAna2020RHCLivetime, bool corrSysts=false, TString options="joint_realData_both_onlyIH", TString plot_type="delta", bool asimov=false)
double dmsq32
const FitSinSqTheta23LowerOctant kFitSinSqTheta23LowerOctant
Definition: FitVars.cxx:17
const double kAna2020RHCLivetime
Definition: Exposures.h:238
void ResetToNominal()
Definition: SystShifts.cxx:141
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
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
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
Float_t e
Definition: plot.C:35
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)
FILE * outfile
Definition: dump_event.C:13
def sign(x)
Definition: canMan.py:197
virtual void SetdCP(const T &dCP)=0
Compare a single data spectrum to the MC + cosmics expectation.
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17