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