joint_fit_2019_contours.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 "Utilities/rootlogon.C"
18 #include "OscLib/IOscCalc.h"
19 
21 //#include "./joint_fit_2019_loader_tools.h"
22 
23 #include "TCanvas.h"
24 #include "TMarker.h"
25 #include "TBox.h"
26 #include "TLatex.h"
27 #include "TColor.h"
28 #include "TGraph.h"
29 #include "TVectorD.h"
30 #include "TF1.h"
31 #include "TLegend.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 DrawHieTag(int hie, bool th13=true, bool vert=false);
45 
46 void FillGraphs(std::vector<TGraph*> g1,std::vector<TGraph*> g2,std::vector<TGraph*> g3,
47  const Int_t kFillColor1, const Int_t kFillColor2, const Int_t kFillColor3,
48  const Int_t kDarkColor, const TString surftype);
49 
50 void DrawReactorConstraint(double mean, double err);
51 
52 // kCMYK analog
54  {
55  Double_t white[9] = { 1, 1, 1, 1, 1, 1, 1, 1, 1 };
56  Double_t Red[9] = { 61./255., 99./255., 136./255., 181./255., 213./255., 225./255., 198./255., 99./255., 61./255.};
57  Double_t Green[9] = { 149./255., 140./255., 96./255., 83./255., 132./255., 178./255., 190./255., 140./255., 149./255};
58  Double_t Blue[9] = { 214./255., 203./255., 168./255., 135./255., 110./255., 100./255., 111./255., 203./255., 214./255.};
59  Double_t Length[9] = { 0.0000, 0.1250, 0.2500, 0.3750, 0.5000, 0.6250, 0.7500, 0.8750, 1.0000};
60 
61  TColor::CreateGradientColorTable(9, Length, Red, Green, Blue, 255, 1.0);
62  gStyle->SetNumberContours(255);
63  }
64 
65 
66 void joint_fit_2019_contours(bool createFile=false,
67  bool corrSysts = true,
68  bool runOnGrid = false,
69  TString options="joint_realData_both_onlyIH",
70  bool dmsqSurf = false,
71  bool th13Surf = false,
72  bool fccorr=false,
73  bool zoomIn = true,
74  bool fillContour = false,
75  bool smooth = true,
76  bool fc_overlay = false)
77 {
78  bool debug = 1;
79  bool nueOnly = options.Contains("nueOnly");
80  bool numuOnly = options.Contains("numuOnly");
81  bool joint = options.Contains("joint");
82  assert (nueOnly || numuOnly || joint);
83 
84  bool FHCOnly = options.Contains("FHCOnly");
85  bool RHCOnly = options.Contains("RHCOnly");
86  bool both = options.Contains("both");
87  assert (FHCOnly || RHCOnly || both);
88 
89  bool fake2018 = options.Contains("fake2018");
90  bool realData = options.Contains("realData");
91 
92  bool onlyNH = options.Contains("onlyNH");
93  bool onlyIH = options.Contains("onlyIH");
94 
95  TString tag;
96  if(FHCOnly) tag="FHCOnly/";
97  if(RHCOnly) tag="RHCOnly/";
98  if(both) tag = "RHC_and_FHC/";
99 
100  auto suffix = options;
101  if(dmsqSurf) suffix+= "_dmsq";
102  if(th13Surf) suffix+= "_th13";
103  if(corrSysts) suffix+="_systs";
104  else suffix+="_stats";
105 
106  TString outdir="/nova/ana/nu_e_ana/Ana2019/Results/"+tag+"contours/FCInputs/";
107  if (runOnGrid) outdir = "./";
108  //outdir = "./";
109  TString outFCfilename (outdir +"contours_FCInput_2019_" + suffix);
110 
111  // outdir="/nova/ana/nu_e_ana/Ana2019/Results/"+tag+"contours/";
112  outdir="./results/blessed_plots/";
113  if(th13Surf) outdir += "th13_delta/";
114  else if (dmsqSurf) outdir += "th23_dmsq/";
115  else outdir += "delta_th23/";
116  outdir += corrSysts?"syst/":"stat/";
117  if (runOnGrid) outdir = "./";
118  //outdir = "./";
119  TString outfilename (outdir + "hist_contours_2019_" + suffix);
120 
121  if(createFile){
122 
123  //////////////////////////////////////////////////
124  // Load Nue and Numu experiments
125  //////////////////////////////////////////////////
126  //need numu only for prestage seeds
127 
128  struct predictions {
129  const string name;
130  const IPrediction * pred;
131  std::pair <TH1D*, double> cos;
132  double pot;
133  };
134 
135  std::vector <predictions> preds;
136  std::vector <predictions> preds_numu_only;
137  std::vector <Spectrum * > data;
138  std::vector <Spectrum * > data_numu_only;
139  std::vector <const IExperiment * > expts;
140  std::vector <const IExperiment * > expts_numu_only;
141 
142  auto calc_fake = DefaultOscCalc();
143  if(fake2018) SetFakeCalc(calc_fake, 0.58, 2.51e-3, 0.17);
144  else if(!realData) {std::cerr << "need setting for data\n"; exit(1);}
145 
146  if(!numuOnly) {
147  if(FHCOnly || both ) {
148  preds.push_back({"nue_fhc", GetNuePrediction2019("combo", DefaultOscCalc(), corrSysts, "fhc", false, runOnGrid), GetNueCosmics2019("fhc", runOnGrid), GetPOT()});
149  }
150  if(RHCOnly || both ) {
151  preds.push_back({"nue_fhc", GetNuePrediction2019("prop", DefaultOscCalc(), corrSysts, "rhc", false, runOnGrid), GetNueCosmics2019("rhc", runOnGrid), GetPOT(false)});
152  }
153 
154  }
155  int nnumu = 4;
156  if(!nueOnly) {
157  if(FHCOnly || both ) {
158  auto numu_preds = GetNumuPredictions2019(nnumu, corrSysts, "fhc", runOnGrid);
159  auto numu_cosmics = GetNumuCosmics2019(nnumu, "fhc", runOnGrid);
160  for (int i = 0; i< nnumu; i++ ){
161  preds.push_back({"numu_fhc_"+std::to_string(i+1), numu_preds[i], numu_cosmics[i], GetPOT()});
162  preds_numu_only.push_back({"numu_fhc_"+std::to_string(i+1), numu_preds[i], numu_cosmics[i], GetPOT()});
163  }
164  }
165  if(RHCOnly || both ) {
166  auto numu_preds = GetNumuPredictions2019(nnumu, corrSysts, "rhc", runOnGrid);
167  auto numu_cosmics = GetNumuCosmics2019(nnumu, "rhc", runOnGrid);
168  for (int i = 0; i< nnumu; i++ ){
169  preds.push_back({"numu_rhc_"+std::to_string(i+1), numu_preds[i], numu_cosmics[i], GetPOT(false)});
170  preds_numu_only.push_back({"numu_rhc_"+std::to_string(i+1), numu_preds[i], numu_cosmics[i], GetPOT(false)});
171  }
172  }
173  }
174 
175  if(realData){
176  if(!numuOnly){
177  if(FHCOnly || both ) data.push_back(GetNueData2019("fhc", runOnGrid));
178  if(RHCOnly || both ) data.push_back(GetNueData2019("rhc", runOnGrid));
179  }
180  if(!nueOnly){
181  if(FHCOnly || both ){
182  auto numu_data = GetNumuData2019(nnumu, "fhc", runOnGrid);
183  data.insert(data.end(),numu_data.begin(), numu_data.end());
184  data_numu_only.insert(data_numu_only.end(),numu_data.begin(), numu_data.end());
185  }
186  if(RHCOnly || both ){
187  auto numu_data = GetNumuData2019(nnumu, "rhc", runOnGrid);
188  data.insert(data.end(),numu_data.begin(), numu_data.end());
189  data_numu_only.insert(data_numu_only.end(),numu_data.begin(), numu_data.end());
190  }
191  }
192  }
193 
194  for(int i = 0; i < int(preds.size()); ++i){
195  double POT = preds[i].pot;
196  auto thisdata = GetFakeData (preds[i].pred, calc_fake, POT, preds[i].cos.first);
197  if(realData) thisdata = data[i];
198  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;
199  expts.push_back(new SingleSampleExperiment(preds[i].pred, *thisdata, preds[i].cos.first, preds[i].cos.second));
200  if(debug){
201  new TCanvas();
202  DataMCComparison (*thisdata, preds[i].pred, DefaultOscCalc());
203  preds[i].cos.first->SetMarkerStyle(34);
204  preds[i].cos.first->SetMarkerColor(kCosmicBackgroundColor);
205  preds[i].cos.first->Draw("hist p same");
206  gPad->Print(outdir + (corrSysts?"syst/":"stat/")+ "debug_predictions_" + suffix + "_" + std::to_string(i) + ".pdf");
207  //gPad->SaveTo();
208  }//end debug
209  }
210 
211  cout<<"Make numu only experiment to get the seeds"<<endl;
212  //make numu only experiment for seeds:
213  for(int i = 0; i < (int) preds_numu_only.size(); ++i){
214  double POT = preds_numu_only[i].pot;
215  auto thisdata = GetFakeData (preds_numu_only[i].pred, calc_fake, POT, preds_numu_only[i].cos.first);
216  if(realData) thisdata = data_numu_only[i];
217  cout<<i<<" "<<preds_numu_only[i].name<<" POT "<<POT<<" tot MC "<<preds_numu_only[i].pred->Predict(calc_fake).Integral(POT)<<" cos "<<preds_numu_only[i].cos.first->Integral()<<" cos er "<<preds_numu_only[i].cos.second<<" analyze data "<<thisdata->Integral(POT)<<endl;
218  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));
219  }
220 
221  //////////////////////////////////////////////////
222  // Add constraints, make experiments
223  //////////////////////////////////////////////////
224  std::cout << "\nCreating multiexperiment\n" << std::endl;
225  if(!th13Surf){
226  std::cout << "Adding WorldReactorConstraint2017\n";
227  expts.push_back(WorldReactorConstraint2017());
228  }
229  if(nueOnly) {
230  std::cout << "Adding Dmsq32ConstraintPDG2017\n";
231  expts.push_back(&kDmsq32ConstraintPDG2017);
232  }
233  std::cout << "Creating Multiexperiment with a total of "
234  << expts.size() << " experiments\n\n" << std::endl;
235  auto exptThis = new MultiExperiment(expts);
236 
237  std::cout << "Creating Multiexperiment of numu only SimpleExp with a total of "
238  << expts_numu_only.size() << " experiments\n\n" << std::endl;
239  auto exptThis_numu_only = new MultiExperiment(expts_numu_only);
240 
241  std::cout<<"Clean verctors"<<std::endl;
242  preds.clear();
243  preds_numu_only.clear();
244 
245  ////////////////////////////////////////////////////////////
246  // Systematics
247  ////////////////////////////////////////////////////////////
248 
249  std::cout << "Systematics for the fit:\n";
250  auto systs = GetJointFitSystematicList(corrSysts, nueOnly, numuOnly, FHCOnly||both, RHCOnly||both, true);
251 
252  std::cout << "\n\nSystematic correlations...\n";
253  if(!nueOnly && ! numuOnly && corrSysts){
254  if(FHCOnly){
255  exptThis->SetSystCorrelations(0, GetCorrelations(true, true));
256  auto notfornumu = GetCorrelations(false, true);
257  for(int i =0; i < nnumu; ++i) exptThis->SetSystCorrelations(i+1, notfornumu);
258  }
259  if(RHCOnly){
260  exptThis->SetSystCorrelations(0, GetCorrelations(true, false));
261  auto notfornumu = GetCorrelations(false, true);
262  for(int i =0; i < nnumu; ++i) exptThis->SetSystCorrelations(i+1, notfornumu);
263  }
264  if(both){
265  exptThis->SetSystCorrelations(0, GetCorrelations(true, true));
266  exptThis->SetSystCorrelations(1, GetCorrelations(true, false));
267  auto notfornumu = GetCorrelations(false, true);
268  for(int i =0; i < 8; ++i) exptThis->SetSystCorrelations(i+2, notfornumu);
269  }
270  }
271  if (nueOnly && corrSysts){
272  if (FHCOnly) exptThis->SetSystCorrelations(0, GetCorrelations(true, true));
273  if (RHCOnly) exptThis->SetSystCorrelations(0, GetCorrelations(true, false));
274  if (both) {
275  exptThis->SetSystCorrelations(0, GetCorrelations(true, true));
276  exptThis->SetSystCorrelations(1, GetCorrelations(true, false));
277  }
278  }
279 
280  std::cout << "Systematics for the numu only fit:\n";
281  auto systs_numu_only = GetJointFitSystematicList(corrSysts, false, true, true, true, true);
282 
283  ////////////////////////////////////////////////////////////
284  // Fit
285  ////////////////////////////////////////////////////////////
286  std::cout << "Starting the fit" << std::endl;
287 
289 
290  SystShifts auxShifts = SystShifts::Nominal();
291  // In case some systs need different seeds
292  std::vector <SystShifts> seedShifts = {};
293 
294  //////////////////////////////////////////////////////////////////////
295  /////////////////////////// Seed controller ////////////////////////////
296  ////////////////////////////////////////////////////////////////////////
297 
298 
299  cout<<"------------------- Start prestage seeds --------------------------"<<endl;
300 
301  double minchi_numu = 1E20;
302  double pre_seed_th23;
303  double pre_seed_dmsq;
304  ResetOscCalcToDefault(calc);
305  auxShifts.ResetToNominal();
306 
307  double maxmix = 0.514; // from the numu results
308  double numu_pre_seedLONH, numu_pre_seedUONH, numu_pre_seedLOIH, numu_pre_seedUOIH, dmsq_numu_pre_seedNH, dmsq_numu_pre_seedIH;
309 
310  for(int hie:{1}){
311  std::cout << "\n\nFinding test best fit " << (hie>0? "NH " : "IH ") << "\n\n";
312  Fitter fitnumu_only(exptThis_numu_only, {&kFitSinSqTheta23, &kFitDmSq32Scaled}, {});
313 
314  auto thisminchi = fitnumu_only.Fit(calc, auxShifts ,
315  {{&kFitDmSq32Scaled,{hie*2.5}},
316  {&kFitSinSqTheta23, {0.4} }}, {},// seedShifts
317  Fitter::kVerbose);
318 
319  if(thisminchi < minchi_numu) minchi_numu = thisminchi;
320  }
321  pre_seed_th23 = kFitSinSqTheta23.GetValue(calc);
322  pre_seed_dmsq = kFitDmSq32Scaled.GetValue(calc);
323 
324  numu_pre_seedLONH = ((pre_seed_th23>maxmix)?(2*maxmix-pre_seed_th23):pre_seed_th23);
325  numu_pre_seedUONH = ((pre_seed_th23>maxmix)?pre_seed_th23:(2*maxmix-pre_seed_th23));
326 
327  ResetOscCalcToDefault(calc);
328  auxShifts.ResetToNominal();
329 
330  cout<<"------------------- End prestage seeds --------------------------"<<endl;
331 
332  struct th23helper{
333  std::vector<double> seeds;
334  const IFitVar * var;
335  TString label;
336  };
337 
338  std::vector <th23helper> th23seeds;
339 
340  th23helper temp;
341  temp.seeds = {0.3};
342  temp.var = kFitSinSqTheta23LowerOctant;
343  temp.label = "LO";
344 
345  //for NH:
346  th23seeds.push_back( { {0.499, numu_pre_seedLONH}, kFitSinSqTheta23LowerOctant, "LO"});
347  th23seeds.push_back( { {0.501, numu_pre_seedUONH}, kFitSinSqTheta23UpperOctant, "UO"});
348  std::vector<double> th23_all_seeds_NH = {numu_pre_seedLONH, 0.5, numu_pre_seedUONH};
349 
350  std::vector<double> delta_seeds = {0, 0.5, 1., 1.5};
351  std::vector<double> th23_old_seeds = {0.45, 0.5, 0.55};
352 
353  /////////////////////////////////////////////////////////////////////////
354  //////////////////////////// Different best fit searches //////////////////
355  ///////////////////////////////////////////////////////////////////////////
356 
357  double minchi23 = 1E20;
358  for(int hie:{-1,1}){
359  for (auto & thseed: th23seeds){
360  std::cout << "\n\nFinding best fit "
361  << thseed.label
362  << (hie>0? " NH " : " IH ")
363  << "\n\n";
364  for (auto const & s:thseed.seeds) std::cout << s << ", ";
365  std::cout << std::endl;
366  cout<<"pre seed from numu is "<<pre_seed_dmsq<<endl;
367  std::vector <const IFitVar*> fitvars = {&kFitDeltaInPiUnits,
368  thseed.var,
371 
372  //if (th13Surf) fitvars = {&kFitDeltaInPiUnits,thseed.var,&kFitDmSq32Scaled};
373 
374  Fitter fit23(exptThis, fitvars, systs);
375  auto thisminchi = fit23.Fit(calc,auxShifts ,{
376  {&kFitDmSq32Scaled,{hie*pre_seed_dmsq}},
377  {thseed.var, thseed.seeds},
378  {&kFitDeltaInPiUnits, delta_seeds}
379  },
380  {}, Fitter::kVerbose
381  );
382  minchi23= min(minchi23, thisminchi);
383  ResetOscCalcToDefault(calc);
384  auxShifts.ResetToNominal();
385  }//end th23
386  }//end hie
387  std::cout << "\nFound overall minchi " << minchi23 << std::endl;
388 
389  cout<<"Check with oldstyle seeds when we didn't split into octant"<<endl;
390  double minchi23test = 1E20;
391  for(int hie:{-1,1}){
392  std::cout << "\n\nFinding test best fit "
393  << (hie>0? "NH " : "IH ")
394  << "\n\n";
395 
396  std::vector <const IFitVar*> fitvars = {&kFitDeltaInPiUnits,
400 
401  //if (th13Surf) fitvars = {&kFitDeltaInPiUnits, &kFitSinSqTheta23, &kFitDmSq32Scaled};
402 
403  Fitter fit23(exptThis, fitvars, systs);
404  auto thisminchi = fit23.Fit(calc,auxShifts ,{
406  {&kFitSinSqTheta23, th23_old_seeds},
407  {&kFitDeltaInPiUnits, delta_seeds}
408  }//, seedShifts
409  );
410  minchi23test= min(minchi23test, thisminchi);
411  ResetOscCalcToDefault(calc);
412  auxShifts.ResetToNominal();
413  }//end hie
414  std::cout << "\nFound overall test minchi " << minchi23test << std::endl;
415 
416  //We'll save to this file, make pretty plots later
417  TFile * outfile = new TFile(outfilename+".root","recreate");
418  //The FC also needs to know about the fitted syst values, save those in
419  //a separate file
420  TFile * outFCfile = new TFile(outFCfilename+".root","recreate");
421 
422  outfile->cd();
423  TVectorD v(1);
424  v[0] = minchi23;
425  cout<<"minchi wrote into the file "<<minchi23<<endl;
426  v.Write("overall_minchi");
427 
428  int steps = 30;//10;// 5;
429  if(dmsqSurf){
430  std::cerr << "\n WARNING Using 20 bins for dmsq surface\n\n";
431  steps = 20;
432  }
433  //Now create all the surfaces
434  for(int hie: {-1, +1}){
435 
436  if((onlyNH && hie<0) || (onlyIH && hie>0)) continue;
437 
438  std::cout << "Starting surface " << (hie>0? "NH ": "IH") << "\n\n";
439  cout<<"the delta seeds are ";
440  for (auto const & s:delta_seeds) std::cout << s << ", ";
441  std::string hieStr = hie>0 ? "NH" : "IH";
442 
443  if(!th13Surf && ! dmsqSurf){
444  ResetOscCalcToDefault(calc);
445  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
446  double low = (nueOnly||RHCOnly)?0:0.3;
447  double high = (nueOnly||RHCOnly)?1:0.7;
448  cout<<"low "<<low<<" high "<<high<<endl;
449  FrequentistSurface surf23(exptThis, calc,
450  &kFitDeltaInPiUnits,steps, 0, 2,
451  &kFitSinSqTheta23, steps, low, high,
453  //{},seedShifts
454  );
455 
456  auto surf1=surf23.ToTH2(minchi23);
457  outfile->cd();
458  surf1->Write((TString)"delta_th23_"+hieStr);
459  surf23.SaveTo(outfile, (TString)"surf_delta_th23_"+hieStr);
460 
461  outFCfile->cd();
462  std::vector<TH2*> profhists = surf23.GetProfiledHists();
463  //Osc parameters are special case
464  profhists[0]->Write((hieStr+"_DmSq32").c_str());
465  profhists[1]->Write((hieStr+"_SinSq2Theta13").c_str());
466  for (int i = 0; i < int(systs.size()); i++){
467  profhists[i+2]->Write((hieStr+"_"+systs[i]->ShortName()).c_str());
468  }
469  }
470  if(!th13Surf && dmsqSurf){
471  calc = DefaultOscCalc();
472  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
473  double low = (nueOnly||RHCOnly)?0.2:0.35;
474  double high = (nueOnly||RHCOnly)?1:0.7;
475  cout<<"low "<<low<<" high "<<high<<endl;
476  FrequentistSurface surf23m(exptThis, calc,
477  &kFitSinSqTheta23, steps,low, high,
478  &kFitDmSq32Scaled, steps,
479  (hie>0?(RHCOnly?2.0:2.2):(RHCOnly?-3.5:-2.8)),(hie>0?(RHCOnly?3.2:2.85):(RHCOnly?-2.0:-2.3)),
481  {{&kFitDeltaInPiUnits, delta_seeds}});
482  auto surf6=surf23m.ToTH2(minchi23);
483  outfile->cd();
484  surf6->Write((TString)"th23_dm32_"+(hie>0?"NH":"IH"));
485  surf23m.SaveTo(outfile, (TString)"surf_th23_dm32_"+(hie>0?"NH":"IH"));
486 
487  outFCfile->cd();
488  std::vector<TH2*> profhists = surf23m.GetProfiledHists();
489  //Osc parameters are special case
490  profhists[0]->Write((hieStr+"_SinSq2Theta13").c_str());
491  profhists[1]->Write((hieStr+"_DeltaCPpi").c_str());
492  for (int i = 0; i < int(systs.size()); i++){
493  profhists[i+2]->Write((hieStr+"_"+systs[i]->ShortName()).c_str());
494  }
495 
496  }
497  if(th13Surf){
498  calc = DefaultOscCalc();
499  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
500 
501  FrequentistSurface surf13(exptThis, calc,
502  &kFitSinSq2Theta13, steps, 0, (RHCOnly?1:0.35),
503  &kFitDeltaInPiUnits,steps, 0, 2,
505  {{&kFitSinSqTheta23, th23_all_seeds_NH}});//, seedShifts);
506  auto surf4 = surf13.ToTH2(minchi23);
507  outfile->cd();
508  surf4->Write((TString)"th13_delta_"+(hie>0?"NH":"IH"));
509  surf13.SaveTo(outfile, (TString)"surf_th13_delta_"+(hie>0?"NH":"IH"));
510 
511  outFCfile->cd();
512  std::vector<TH2*> profhists = surf13.GetProfiledHists();
513  //Osc parameters are special case
514  profhists[0]->Write((hieStr+"_SinSqTheta23").c_str());
515  profhists[1]->Write((hieStr+"_DmSq32").c_str());
516  for (int i = 0; i < int(systs.size()); i++){
517  profhists[i+2]->Write((hieStr+"_"+systs[i]->ShortName()).c_str());
518  }
519 
520  }
521  } //end hie
522 
523  delete outfile;
524  std::cout << "\nSurfaces saved to " << outfilename << std::endl;
525  return;
526  }//end createfile
527 
528 //////////////////////////////////////////////////////////////////////
529 // Make plots
530 //////////////////////////////////////////////////////////////////////
531  else{ //plot only
532 
533  TString infilename = "/nova/ana/nu_e_ana/Ana2019/Results/RHC_and_FHC/contours/";
534  if(th13Surf) infilename += "th13_delta/";
535  else if (dmsqSurf) infilename += "th23_dmsq/";
536  else infilename += "delta_th23/";
537  infilename += corrSysts?"syst/":"stat/";
538  infilename += "hist_contours_2019_";
539  infilename += options;
540  if (dmsqSurf) infilename += "_dmsq";
541  infilename += corrSysts?"_systs.root":"_stats.root";
542 
543  TFile * infile = new TFile (infilename,"read");
544  TFile * official =0;
545  if(corrSysts && options.Contains("joint") && fccorr && 0) {
546  official = new TFile("NOvA_nue_official_results.root","recreate");
547  TH2* axes = new TH2F("axes",";#delta_{CP}/#pi;sin^{2}#theta_{23}",100,0,2,100,0.225,0.725);
548  //CenterTitles(axes);
549  axes->Write();
550  }
551 
552  auto mins =* (TVectorD*)infile->Get("overall_minchi");
553  double minchi23 = mins[0];
554 
555  std::vector <TString> surfNames;// = {"delta_th23_"};<
556  if(!th13Surf && ! dmsqSurf)surfNames = {"delta_th23_"};
557  if (dmsqSurf) surfNames.push_back("th23_dm32_");
558  if (th13Surf) surfNames.push_back("th13_delta_");
559  for(TString surfName:surfNames){
560 
561  for (int hie:{-1,+1}){
562 
563  if((onlyNH && hie<0) || (onlyIH && hie>0)) continue;
564 
565  //Setup canvas
566  if(surfName.Contains("delta_th23")) {
567  if (zoomIn) DrawContourCanvas(surfName, 0., 2., (RHCOnly?0.15:0.275), (RHCOnly?0.8:0.725));
568  else DrawContourCanvas(surfName, 0., 2., 0., 1.);
569  }
570  if(surfName.Contains("th23_dm32")) {
571  if (hie>0) DrawContourCanvas(surfName, ((nueOnly||RHCOnly)?0.2:0.3), ((nueOnly||RHCOnly)?0.8:0.7), (RHCOnly?1.8:2.05), (RHCOnly?3.3:2.85));
572  else DrawContourCanvas(surfName, ((nueOnly||RHCOnly)?0.2:0.3), ((nueOnly||RHCOnly)?0.8:0.7), (RHCOnly?-3.5:-2.95), (RHCOnly?-2.0:-2.25));
573  }
574  if(surfName.Contains("th13_delta")) {
575  if (zoomIn) DrawContourCanvas(surfName, 0., 0.38, 0, 2);
576  else DrawContourCanvas(surfName, 0., 1.0, 0, 2);
577  }
578  //Load the surface; fill and FC corrections if available
579  auto surf = *FrequentistSurface::LoadFrom(infile, "surf_"+surfName+(hie>0?"NH":"IH"));
580  TH2 * hsurf = (TH2*)infile->Get(surfName+(hie>0?"NH":"IH"));
581  TH2 * surf_1Sigma = Gaussian68Percent2D(surf);
582  TH2 * surf_2Sigma = Gaussian2Sigma2D(surf);
583  TH2 * surf_3Sigma = Gaussian3Sigma2D(surf);
584  TH2 * fc_surf_1Sigma = new TH2F();
585  TH2 * fc_surf_2Sigma = new TH2F();
586  TH2 * fc_surf_3Sigma = new TH2F();
587  if(fccorr){
588  // if(!joint || !corrSysts || decomp != "combo" || dmsqSurf || th13Surf)
589  // { std::cerr << "\nWe don't have FC\n\n"; exit(1);}
590 
591  TString fcFolder="/nova/ana/nu_e_ana/Ana2019/FC/";
592  TString fcFileName;
593  if(dmsqSurf) fcFileName = "contours/ssth23dmsq/ssth23dmsq_";
594  else if(!dmsqSurf && !th13Surf)
595  fcFileName = "contours/deltassth23/deltassth23_";
596  fcFileName += hie > 0? "nh":"ih";
597  fcFileName += ".root";
598 
599  FCSurface *fcXH = 0;
600  auto fccol = FCCollection::LoadFromFile((fcFolder + fcFileName).Data()).release();
601  int NX = hsurf->GetXaxis()->GetNbins();
602  int NY = hsurf->GetYaxis()->GetNbins();
603  double minX = hsurf->GetXaxis()->GetBinLowEdge(1);
604  double maxX = hsurf->GetXaxis()->GetBinUpEdge(NX);
605  double minY = hsurf->GetYaxis()->GetBinLowEdge(1);
606  double maxY = hsurf->GetYaxis()->GetBinUpEdge(NY);
607 
608  std::cout << "Making FCSurface: \n";
609  std::cout << "-------------------------------\n";
610  std::cout << "NX: " << NX;
611  std::cout << "\tNY: " << NY;
612  std::cout << "\tminX: " << minX;
613  std::cout << "\tmaxX: " << maxX;
614  std::cout << "\tminY: " << minY;
615  std::cout << "\tmaxY: " << maxY;
616  std::cout << "\n------------------------------\n";
617 
618  fcXH = new FCSurface(NX, minX, maxX, NY, minY, maxY);
619  std::string plot_name = "deltassth23";
620  if(dmsqSurf) plot_name = "ssth23dmsq32";
621  if(th13Surf) plot_name = "deltath13";
622 
623  fcXH->Add(*fccol, plot_name);
624  fc_surf_1Sigma = fcXH->SurfaceForSignificance(0.6827);
625  fc_surf_2Sigma = fcXH->SurfaceForSignificance(0.9545);
626  fc_surf_3Sigma = fcXH->SurfaceForSignificance(0.9973);
627  if(smooth){
628  fc_surf_1Sigma->Smooth();
629  fc_surf_2Sigma->Smooth();
630  if(dmsqSurf || th13Surf)
631  fc_surf_3Sigma->Smooth();
632  else
633  CylindricalSmoothing(fc_surf_3Sigma);
634  }
635  }
636  Color_t k2SigmaConfidenceColor = ((hie>0)? cnh:cih);
637  Int_t kFillColor1 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.99);
638  Int_t kFillColor2 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.5);
639  Int_t kFillColor3 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.2);
640  Int_t kFCFillColor1 = TColor::GetColorTransparent(k2SigmaConfidenceColor, fc_overlay? 0.65: 0.75);
641  Int_t kFCFillColor2 = TColor::GetColorTransparent(k2SigmaConfidenceColor, fc_overlay? 0.4: 0.38);
642  Int_t kFCFillColor3 = TColor::GetColorTransparent(k2SigmaConfidenceColor, fc_overlay? 0.18: 0.18);
643  Int_t kDarkColor = kGray+3;
644 
645  if(fillContour){
646  if(!fccorr){
647  auto g1 = surf.GetGraphs(surf_1Sigma,minchi23);
648  auto g2 = surf.GetGraphs(surf_2Sigma,minchi23);
649  auto g3 = surf.GetGraphs(surf_3Sigma,minchi23);
650  FillGraphs(g1, g2, g3,
651  kFillColor1, kFillColor2, kFillColor3,
652  kDarkColor, suffix + surfName + (hie>0?"NH":"IH"));
653  }
654  else {
655  auto fc_g1 = surf.GetGraphs(fc_surf_1Sigma,minchi23);
656  auto fc_g2 = surf.GetGraphs(fc_surf_2Sigma,minchi23);
657  auto fc_g3 = surf.GetGraphs(fc_surf_3Sigma,minchi23);
658  FillGraphs(fc_g1, fc_g2, fc_g3,
659  kFillColor1, kFillColor2, kFillColor3,
660  kDarkColor, suffix + surfName + (hie>0?"NH":"IH")+ "fccorr");
661  }
662  }
663  else{
664  cout<<"minchi from file is "<<minchi23<<endl;
665  if(!fccorr || fc_overlay){
666 
667  if(fc_overlay) {
668  auto g1 = surf.GetGraphs(surf_1Sigma,minchi23);
669  auto g2 = surf.GetGraphs(surf_2Sigma,minchi23);
670  auto g3 = surf.GetGraphs(surf_3Sigma,minchi23);
671  FillGraphs(g1, g2, g3,
672  kFillColor1, kFillColor2, kFillColor3,
673  kDarkColor, suffix + surfName + (hie>0?"NH":"IH"));
674  }
675  surf.DrawContour(surf_1Sigma, kSolid, kFillColor1,minchi23);
676  surf.DrawContour(surf_2Sigma, k2SigmaConfidenceStyle, kFillColor2,minchi23);
677  surf.DrawContour(surf_3Sigma, kSolid, kFillColor3,minchi23);
678  }
679  if(fccorr) {
680  surf.DrawContour(fc_surf_1Sigma, kSolid, kFCFillColor1, minchi23);
681  surf.DrawContour(fc_surf_2Sigma, k2SigmaConfidenceStyle, kFCFillColor2, minchi23);
682  surf.DrawContour(fc_surf_3Sigma, kSolid, kFCFillColor3, minchi23);
683  }
684  }
685 
686  gPad->RedrawAxis();
687 
688  bool bestfit = 1;
689  double delta, th23, dmsq32, th13;
690  if(bestfit && both && corrSysts && hie>0){
691  if(!th13Surf){
692  delta = 1.999; /// Insert new values or make them read from the file
693  th23 = 0.5649;
694  dmsq32 = 2.483;
695  th13 = 0.08189;
696  }
697  else{delta = 0.678516; th13 = 0.0957894;}
698  TMarker bf;
699  bf.SetMarkerColor(kDarkColor);
700  // bf.SetMarkerStyle(kFullStar);
701  bf.SetMarkerStyle(43);
702  bf.SetMarkerSize(2);
703  if(!dmsqSurf && ! th13Surf){
704  bf.SetX(delta);
705  bf.SetY(th23);
706  bf.DrawClone();
707  }
708  if(dmsqSurf){
709  bf.SetX(th23);
710  bf.SetY(dmsq32);
711  bf.DrawClone();
712  }
713  if(th13Surf){
714  bf.SetX(th13);
715  bf.SetY(delta);
716  bf.DrawClone();
717  }
718  }
719  auto leg = ContourLegend(hie, fillContour, fccorr,
720  kFillColor1,kFillColor2,kFillColor3,kDarkColor,
721  bestfit && both);
722 
723  if(th13Surf) {
724  DrawReactorConstraint(0.082,0.004);
725  leg->SetFillStyle(1001);
726  leg->SetNColumns(1);
727  leg->SetX1(0.7); leg->SetX2(0.8);
728  leg->SetY1(0.4); leg->SetY2(0.85);
729  if(!fccorr) leg->SetHeader("No FC");
730 
731  TGraph * dummy = new TGraph;
732  dummy->SetLineWidth(2);
733  dummy->SetLineColor(kGray+2);
734  dummy->SetFillColor(kGray);
735  leg->AddEntry(dummy,"#splitline{Reactor}{68%C.L.}","f");
736  }
737  leg->Draw();
738  Preliminary();
739  DrawHieTag(hie,false);
740 
741 
742  for(TString ext: {".png", ".pdf",".root"}){
743  gPad->Print(outdir + "contour_" + suffix + "_" +
744  surfName +
745  (hie>0?"NH":"IH") + (fccorr?"_fccorr":"") +
746  ((fccorr && smooth)?"_smooth":"") +
747  (zoomIn?"_zoom":"") +
748  (fc_overlay?"_overlay":"") +
749  ext);
750  }//end ext
751 
752  debug = true; // I don't want 10 million plots
753  if(debug && !fccorr){
754  auto hists = surf.GetProfiledHists();
755  int nSysts = corrSysts ? 51 : 0;
756  for (int i = 0; i < nSysts; ++i){
757  hists.push_back((TH2*) infile->Get((TString) "surf_" + surfName +
758  (hie > 0 ? "NH" : "IH") +
759  "/profHists/hist" + std::to_string(i)));
760  }
761  hists.push_back(surf.ToTH2(minchi23));
762 
763  for(uint i = 0; i < hists.size(); ++i){
764  if(!hists[i]) continue;
765  if(TString(hists[i]->GetTitle()).Contains("delta")) {
767  hists[i]->GetZaxis()->SetRangeUser(0,2);
768  }
769  else gStyle->SetPalette(87);
770  if( i > 3) {
771  double minz = hists[i]->GetBinContent(hists[i]->GetMinimumBin());
772  double maxz = hists[i]->GetBinContent(hists[i]->GetMaximumBin());
773  double lim = std::max(fabs(minz), fabs(maxz));
774  hists[i]->GetZaxis()->SetRangeUser(-lim,+lim);
775  }
776  if(i==(hists.size()-1)) {
777  gStyle->SetPalette(56);
778  hists[i]->GetZaxis()->SetRangeUser(0.,11.83);
779  }
780 
781  hists[i]->Draw("colz");
782  hists[i]->GetYaxis()->SetTitleOffset(0.8);
783  gPad->SetRightMargin(0.14);
784  surf.DrawContour(surf_1Sigma, kSolid, kDarkColor,minchi23);
785  surf.DrawContour(surf_2Sigma, kSolid, kDarkColor,minchi23);
786  surf.DrawContour(surf_3Sigma, kSolid, kDarkColor,minchi23);
787  gPad->Print(outdir + "debug_contour_" + suffix + "_" +
788  surfName +
789  (hie > 0 ? "NH":"IH") + (fccorr?"_fccorr":"") +
790  (zoomIn ? "_zoom" : "") +
791  (std::to_string(i)) + ".pdf");
792 
793  }//end profiles
794  }//end debug
795  }//end hie
796  }//end surf
797  } // end plot only
798 }//end all
799 
800 
801 void DrawHieTag(int hie,bool th13, bool vert){
802  if(hie>0){
803  TLatex * ltxNH = new TLatex(0.16,0.89,"#bf{NH}");
804  ltxNH->SetTextColor(cnh);
805  ltxNH->SetNDC();
806  ltxNH->SetTextAlign(22);
807  ltxNH->SetTextSize(kBlessedLabelFontSize);
808  if(th13) ltxNH->Draw();
809  else if(!vert) ltxNH->DrawLatex(.85,.20,"#bf{NH}");
810  else ltxNH->DrawLatex(.85,.85,"#bf{NH}");
811  }
812  else{
813  TLatex * ltxIH = new TLatex (0.16,0.46,"#bf{IH}");
814  ltxIH->SetTextColor(cih_dark);
815  ltxIH->SetNDC();
816  ltxIH->SetTextAlign(22);
817  ltxIH->SetTextSize(kBlessedLabelFontSize);
818  if(th13) ltxIH->Draw();
819  else if(!vert) ltxIH->DrawLatex(.85,.20,"#bf{IH}");
820  else ltxIH->DrawLatex(.85,.85,"#bf{IH}");
821  }
822 }
823 
824 
825 void FillGraphs(std::vector<TGraph*> g1,std::vector<TGraph*> g2,std::vector<TGraph*> g3,
826  const Int_t kFillColor1, const Int_t kFillColor2, const Int_t kFillColor3,
827  const Int_t kDarkColor, const TString surftype)
828 {
829  bool isJoint = surftype.Contains("joint");
830  bool isBoth = surftype.Contains("both");
831  bool isRHC = surftype.Contains("RHCOnly");
832  bool isNH = surftype.Contains("NH");
833  bool isSysts = surftype.Contains("systs");
834  bool isFccorr = surftype.Contains("fccorr");
835 
836  if (surftype.Contains("delta_th23")){
837  cout<<"delta"<<endl;
838  if(isJoint && isNH ){
839  cout<<"joint nh"<<endl;
840  if(isRHC){
841  JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f");
842  JoinGraphs(g2[0], g2[1], kWhite)->Draw("f");
843  JoinGraphs(g2[1], g2[2], kWhite)->Draw("f");
844  JoinGraphs(g2[0], g2[1], kFillColor2)->Draw("f same");
845  JoinGraphs(g2[1], g2[2], kFillColor2)->Draw("f same");
846  }
847  else{
848  JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f");
849  JoinGraphs(g3[0], g3[1], kWhite)->Draw("l");
850  JoinGraphs(g2[0], g2[1], kWhite)->Draw("f");
851  JoinGraphs(g2[0], g2[1], kFillColor2)->Draw("f");
852  JoinGraphs(g2[0], g2[1], kWhite)->Draw("l");
853  JoinGraphs(g1[0], g1[1], kFillColor1)->Draw("c f");
854  JoinGraphs(g1[0], g1[1], kWhite)->Draw("c");
855  }
856  // if(isBoth){
857  // for (auto & g:g1) { g->SetFillColor(kWhite); g->DrawClone("c f");}
858  // for (auto & g:g1) { g->SetFillColor(kFillColor1); g->Draw("c f");}
859  // for (auto & g:g1) { g->Draw("c");}
860  // }
861  }
862 
863  if(isJoint && !isNH){
864  if(isRHC) JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f");
865  else{
866  //JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f c");
867  g3[1]->SetFillColor(kFillColor3); g3[1]->Draw("f c");
868  g3[0]->SetFillColor(kFillColor3); g3[0]->Draw("f c");
869  for (auto &g:g3) {
870  g->SetLineWidth(3);
871  g->Draw("c");
872  }
873  }
874  for (auto &g:g2) {
875  g->SetFillColor(kWhite); g->DrawClone("f c");
876  g->SetFillColor(kFillColor2); g->Draw("f c");
877  g->SetLineColor(kDarkColor);
878  g->SetLineWidth(3);
879  g->Draw("c");
880  }
881  for (auto &g:g1) {
882  g->SetFillColor(kWhite); g->DrawClone("f c");
883  g->SetFillColor(kFillColor1); g->Draw("f c");
884  g->SetLineColor(kDarkColor);
885  g->SetLineWidth(3);
886  g->Draw("c");
887  }
888  for(auto & gs:{g3, g2, g1}){
889  for(auto & g:gs) {
890  g->SetLineColor(kDarkColor);
891  g->SetLineWidth(3);
892  g->Draw("c");
893  }
894  }
895  }
896  }
897 
898  if (surftype.Contains("th23_dm32")){
899  for (auto &g:g3) {
900  g->SetFillColor(kFillColor3); g->Draw("cf");
901  }
902  for (auto &g:g2) {
903  g->SetFillColor(kWhite); g->DrawClone("cf");
904  g->SetFillColor(kFillColor2); g->Draw("cf");
905  }
906  for (auto &g:g1) {
907  g->SetFillColor(kWhite); g->DrawClone("cf");
908  g->SetFillColor(kFillColor1); g->Draw("cf");
909  }
910  for (auto & gs:{g3,g2,g1}){
911  for (auto & g:gs) {
912  g->SetLineColor(kDarkColor);
913  g->SetLineWidth(3);
914  g->Draw("c");
915  }
916  }
917 
918  }
919 
920  if (surftype.Contains("th13_delta")) {
921  if(isBoth) {
922  JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f");
923  if(isNH) {
924  JoinGraphs(g2[0], g2[1], kWhite)->Draw("f");
925  JoinGraphs(g2[0], g2[1], kFillColor2)->Draw("f");
926  }
927  else{
928  for (auto & g:g2) { g->SetFillColor(kWhite); g->DrawClone("f");}
929  for (auto & g:g2) { g->SetFillColor(kFillColor2); g->Draw("f");}
930  }
931  for (auto & g:g1) { g->SetFillColor(kWhite); g->DrawClone("f");}
932  for (auto & g:g1) { g->SetFillColor(kFillColor1); g->Draw("f");}
933  }
934  else{
935  //for (auto &g:g3) {
936  // g->SetFillColor(kFillColor3); g->Draw("f");}
937  //JoinGraphs(g3[0], g3[1], kWhite)->Draw("f");
938  JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f c");
939  JoinGraphs(g2[0], g2[1], kWhite)->Draw("f c");
940  JoinGraphs(g2[0], g2[1], kFillColor2)->Draw("f c");
941  /*for (auto &g:g2) {
942  g->SetFillColor(kWhite); g->DrawClone("f");
943  g->SetFillColor(kFillColor2); g->Draw("f");
944  }*/
945  JoinGraphs(g1[0], g1[1], kWhite)->Draw("f c");
946  JoinGraphs(g1[0], g1[1], kFillColor1)->Draw("f c");
947  }
948 
949  for (auto & gs:{g3,g2,g1}){
950  for (auto & g:gs) {
951  g->SetLineColor(kDarkColor);
952  g->SetLineWidth(3);
953  g->Draw("l");
954  }
955  }
956  }
957 
958 }
959 
960 void DrawReactorConstraint(double mean, double err)
961 {
962  TBox* box = new TBox(mean-err, 0, mean+err, 2);
963  box->SetFillColorAlpha(kGray,0.7);
964 
965  TLine* linLo = new TLine(mean-err, 0, mean-err, 2);
966  linLo->SetLineColor(kGray+2);
967  linLo->SetLineWidth(2);
968 
969  TLine* linHi = new TLine(mean+err, 0, mean+err, 2);
970  linHi->SetLineColor(kGray+2);
971  linHi->SetLineWidth(2);
972 
973  box->Draw();
974  linLo->Draw();
975  linHi->Draw();
976  gPad->RedrawAxis();
977 }
978 
979 // ssth23dmsq32 NH
980 // cafe joint_fit_2019_contours.C false true false joint_realData_both_onlyNH true false true true false true false
981 // ssth23dmsq32 IH
982 // cafe joint_fit_2019_contours.C false true false joint_realData_both_onlyIH true false true true false true false
983 
984 // deltassth23 NH
985 // cafe joint_fit_2019_contours.C false true false joint_realData_both_onlyNH false false true true false true false
986 // deltassth23 IH
987 // cafe joint_fit_2019_contours.C false true false joint_realData_both_onlyIH false false true true false true false
Spectrum * GetFakeData(const IPrediction *pred, osc::IOscCalc *calc, const double pot, TH1D *cosmics=0)
T max(const caf::Proxy< T > &a, T b)
void FillGraphs(std::vector< TGraph * > g1, std::vector< TGraph * > g2, std::vector< TGraph * > g3, const Int_t kFillColor1, const Int_t kFillColor2, const Int_t kFillColor3, const Int_t kDarkColor, const TString surftype)
double minY
Definition: plot_hist.C:9
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
void SetPaletteDeltaCyclicNew()
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 Color_t kCosmicBackgroundColor
Definition: Style.h:41
double th23
Definition: runWimpSim.h:98
double delta
Definition: runWimpSim.h:98
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
TH2 * Gaussian68Percent2D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 2D in gaussian approximation.
OStream cerr
Definition: OStream.cxx:7
double maxY
Definition: plot_hist.C:10
TH2 * Gaussian2Sigma2D(const FrequentistSurface &s)
Up-value surface for 2 sigma confidence in 2D in gaussian approximation.
void CylindricalSmoothing(TH2 *h)
TString hists[nhists]
Definition: bdt_com.C:3
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)
virtual void SetDmsq32(const T &dmsq32)=0
const char * label
tuple white
Definition: rootlogon.py:71
const XML_Char const XML_Char * data
Definition: expat.h:268
Log-likelihood scan across two parameters.
const XML_Char * s
Definition: expat.h:262
string outfilename
knobs that need extra care
string infile
void DrawReactorConstraint(double mean, double err)
void DrawHieTag(int hie, bool th13=true, bool vert=false)
TH2 * Gaussian3Sigma2D(const FrequentistSurface &s)
Up-value surface for 3 sigma confidence in 2D in gaussian approximation.
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
virtual T GetDmsq32() const
Definition: IOscCalc.h:91
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
void Add(const FCPoint &pt, std::string plot)
Definition: FCSurface.cxx:39
void SetFakeCalc(osc::IOscCalcAdjustable *calc, double ssth23=-5, double dmsq32=-5, double dCP_Pi=-5)
void Preliminary()
void DrawContourCanvas(TString surfName, double minx=0, double maxx=2, double miny=0, double maxy=1)
TGraph * JoinGraphs(TGraph *a, TGraph *b, int fillcol)
Join graphs and set a fill color. Useful for contours.
Definition: Plots.cxx:1529
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
Spectrum * GetNueData2019(std::string beam, bool GetFromUPS=false)
TH2 * ToTH2(double minchi=-1) const
Definition: ISurface.cxx:161
TLegend * ContourLegend(int hie, bool fillContour, bool fccorr, Int_t kFillColor1, Int_t kFillColor2, Int_t kFillColor3, Int_t kDarkColor, bool bestFit)
Combine multiple component experiments.
const Style_t k2SigmaConfidenceStyle
Definition: Style.h:71
std::vector< std::pair< TH1D *, double > > GetNumuCosmics2019(const int nq=4, std::string beam="fhc", bool GetFromUPS=false, bool NERSC=false)
double dmsq32
std::vector< Spectrum * > GetNumuData2019(const int nq=4, std::string beam="fhc", bool GetFromUPS=false)
const FitSinSqTheta23LowerOctant kFitSinSqTheta23LowerOctant
Definition: FitVars.cxx:17
void ResetToNominal()
Definition: SystShifts.cxx:144
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
exit(0)
void joint_fit_2019_contours(bool createFile=false, bool corrSysts=true, bool runOnGrid=false, TString options="joint_realData_both_onlyIH", bool dmsqSurf=false, bool th13Surf=false, bool fccorr=false, bool zoomIn=true, bool fillContour=false, bool smooth=true, bool fc_overlay=false)
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
const Float_t kBlessedLabelFontSize
Definition: Style.h:90
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
TH2 * SurfaceForSignificance(double sig)
Definition: FCSurface.cxx:103
Float_t e
Definition: plot.C:35
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
static std::unique_ptr< FrequentistSurface > LoadFrom(TDirectory *dir, const std::string &name)
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)
double th13
Definition: runWimpSim.h:98
static std::unique_ptr< FCCollection > LoadFromFile(const std::string &wildcard)
const std::string outdir
FILE * outfile
Definition: dump_event.C:13
Pseudo-experiments, binned to match a log-likelihood surface.
Definition: FCSurface.h:14
std::pair< TH1D *, double > GetNueCosmics2019(std::string beam, bool GetFromUPS=false, bool NERSC=false)
Compare a single data spectrum to the MC + cosmics expectation.
surf
Definition: demo5.py:35
enum BeamMode string
unsigned int uint