joint_fit_2017_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 "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 
39 using namespace ana;
40 
41 void DrawHieTag(int hie, bool th13=true, bool vert=false);
42 
43 TLegend * ContourLegend(int hie, bool fillContour, bool fccorr,
44  Int_t kFillColor1, Int_t kFillColor2,
45  Int_t kFillColor3, Int_t kDarkColor,
46  bool bestFit);
47 void FillGraphs(std::vector<TGraph*> g1,std::vector<TGraph*> g2,std::vector<TGraph*> g3,
48  const Int_t kFillColor1, const Int_t kFillColor2, const Int_t kFillColor3,
49  const Int_t kDarkColor, const TString surftype);
50 
51 void DrawContourCanvas (TString surfName, double minx = 0,
52  double maxx = 2, double miny = 0, double maxy = 1){
53  auto c1 = new TCanvas(ana::UniqueName().c_str());
54  c1->SetFillStyle(0);
55  c1->SetLeftMargin(0.14);
56  c1->SetBottomMargin(0.15);
57  TH2* axes = new TH2F();
58  TString title;
59  if(surfName.Contains("delta_th23")) title=";#delta_{CP};sin^{2}#theta_{23}";
60  if(surfName.Contains("th13_delta")) title=";sin^{2}2#theta_{13};#delta_{CP}/#pi";
61  if(surfName.Contains("th23_dm32")) title=";sin^{2}#theta_{23};#Deltam^{2}_{32} (10^{-3}eV^{2})";
62  axes = new TH2F("",title,100,minx,maxx,100,miny,maxy);
63 
64  CenterTitles(axes);
65  axes->Draw();
66  axes->GetXaxis()->SetTitleSize(kBlessedTitleFontSize);
67  axes->GetYaxis()->SetTitleSize(kBlessedTitleFontSize);
68  axes->GetXaxis()->SetLabelSize(kBlessedLabelFontSize);
69  axes->GetYaxis()->SetLabelSize(kBlessedLabelFontSize);
70  axes->GetXaxis()->SetTitleOffset(0.8);
71  axes->GetYaxis()->SetTitleOffset(0.8);
72  TGaxis::SetMaxDigits(3);
73  if(surfName.Contains("th23_dm32")) {
74  axes->GetYaxis()->SetDecimals() ;
75  axes->GetYaxis()->SetTitleOffset(0.85);
76  axes->GetYaxis()->SetLabelOffset(0.002);
77  axes->GetYaxis()->SetLabelSize(0.058);
78  axes->GetYaxis()->SetTitleSize(0.078);
79  }
80  if(surfName.Contains("delta_th23")) XAxisDeltaCPLabels(axes);
81  c1->RedrawAxis();
82 }
83 
85  bool createFile=false,
86  bool corrSysts = false,
87  TString options="joint_fakeNHLO",
88  bool dmsqSurf = false,
89  bool th13Surf = false,
90  bool fccorr=false,
91  bool zoomIn = true,
92  bool fillContour = false,
93  bool smooth = true)
94 {
95  bool debug = 1;
96  bool nueOnly = options.Contains("nueOnly");
97  bool numuOnly = options.Contains("numuOnly");
98  bool joint = options.Contains("joint");
99  assert (nueOnly || numuOnly || joint);
100 
101  bool fakeNHLO = options.Contains("fakeNHLO");
102  bool fakeNHUO = options.Contains("fakeNHUO");
103  bool fake2017 = options.Contains("fake2017");
104  bool realData = options.Contains("realData");
105 
106  bool onlyNH = options.Contains("onlyNH");
107  bool onlyIH = options.Contains("onlyIH");
108 
109  auto suffix = options;
110  suffix += decomp;
111  if(dmsqSurf) suffix+= "_dmsq";
112  if(th13Surf) suffix+= "_th13";
113  if(corrSysts) suffix+="_systs";
114  else suffix+="_stats";
115 
116  TString outdir="/nova/ana/nu_e_ana/Ana2017/Results/contours/";
117  if(th13Surf) outdir += "th13_delta/";
118  else if (dmsqSurf) outdir += "th23_dmsq/";
119  else outdir += "delta_th23/";
120 
121  //TString outdir = "contours/";
122 
123  TString outfilename (outdir + "hist_contours_2017_" + suffix);
124  TString outFCfilename (outdir + "contours_FCInput_2017_" + decomp + ".root");
125 
126  if(createFile){
127 
128  //////////////////////////////////////////////////
129  // Load Nue and Numu experiments
130  //////////////////////////////////////////////////
131 
132  std::vector <const IPrediction * > preds;
133  std::vector <std::pair <TH1D*, double > > cosmics;
134  std::vector <Spectrum * > data;
135  std::vector <const IExperiment * > expts;
136 
137  auto calc_fake = DefaultOscCalc();
138  if(fakeNHLO) SetFakeCalc(calc_fake, 0.404, 2.7e-3, 1.48);
139  else if(fakeNHUO) SetFakeCalc(calc_fake, 0.623, 2.7e-3, 0.74);
140  else if(fake2017) SetFakeCalc(calc_fake, 0.56, 2.43e-3, 1.21);
141  else if(!realData) {std::cerr << "need setting for data\n"; exit(1);}
142 
143  if(!numuOnly) {
144  preds.push_back(GetNuePrediction2017(decomp, DefaultOscCalc(), corrSysts));
145  cosmics.push_back(GetNueCosmics2017());
146  }
147  int nnumu = 4;
148  if(!nueOnly) {
149  auto numu_preds = GetNumuPredictions2017(nnumu, corrSysts);
150  preds.insert(preds.end(),numu_preds.begin(), numu_preds.end());
151  auto numu_cosmics = GetNumuCosmics2017(nnumu);
152  cosmics.insert(cosmics.end(),numu_cosmics.begin(), numu_cosmics.end());
153  }
154  if(realData){
155  if(!numuOnly){
156  data.push_back(GetNueData2017());
157  }
158  if(!nueOnly){
159  auto numu_data = GetNumuData2017(nnumu);
160  data.insert(data.end(),numu_data.begin(), numu_data.end());
161  }
162  }
163  for(int i = 0; i < int(preds.size()); ++i){
164  auto thisdata = GetFakeData (preds[i],calc_fake, kAna2017POT,
165  cosmics[i].first);
166  if(realData) thisdata = data[i];
167  expts.push_back(new SingleSampleExperiment(preds[i],*thisdata,
168  cosmics[i].first,cosmics[i].second));
169 
170  if(debug){
171  new TCanvas();
172  DataMCComparison (*thisdata, preds[i], DefaultOscCalc());
173  cosmics[i].first->SetMarkerStyle(34);
174  cosmics[i].first->SetMarkerColor(kCosmicBackgroundColor);
175  cosmics[i].first->Draw("hist p same");
176  gPad->Print(outdir + "debug_predictions_" + suffix + "_" + std::to_string(i) + ".pdf");
177  }//end debug
178  }
179 
180  //////////////////////////////////////////////////
181  // Add constraints, make experiments
182  //////////////////////////////////////////////////
183  std::cout << "\nCreating multiexperiment\n" << std::endl;
184  if(!th13Surf){
185  std::cout << "Adding WorldReactorConstraint2017\n";
186  expts.push_back(WorldReactorConstraint2017());
187  }
188  if(nueOnly) {
189  std::cout << "Adding Dmsq32ConstraintPDG2017\n";
190  expts.push_back(&kDmsq32ConstraintPDG2017);
191  }
192  std::cout << "Creating Multiexperiment with a total of "
193  << expts.size() << " experiments\n\n" << std::endl;
194  auto exptThis = new MultiExperiment(expts);
195 
196  ////////////////////////////////////////////////////////////
197  // Systematics
198  ////////////////////////////////////////////////////////////
199  std::cout << "Systematics for the fit:\n";
200  auto systs = GetJointFitSystematicList(corrSysts, !numuOnly, !nueOnly, true);
201 
202 // std::cout << "HACK MECQ0\n";
203 // systs = {kMECq0ShapeSyst2017};
204 
205  std::cout << "\n\nSystematic correlations...\n";
206  if(!nueOnly && ! numuOnly && corrSysts){
207  exptThis->SetSystCorrelations(0, GetCorrelations(true));
208  auto notfornumu = GetCorrelations(false);
209  for(int i =0; i < nnumu; ++i) exptThis->SetSystCorrelations(i+1, notfornumu);
210  }
211 
212  ////////////////////////////////////////////////////////////
213  // Fit
214  ////////////////////////////////////////////////////////////
215  std::cout << "Starting the fit" << std::endl;
216 
218 
219  std::vector <const IFitVar *> fitvars = {&kFitDeltaInPiUnits,&kFitSinSqTheta23,
221  if (th13Surf) fitvars = {&kFitDeltaInPiUnits,&kFitSinSqTheta23,&kFitDmSq32Scaled};
222 
223  //Find the best fit
224  MinuitFitter fit23(exptThis,fitvars,systs);
225 
226  SystShifts auxShifts = SystShifts::Nominal();
227  // In case some systs need different seeds
228  std::vector <SystShifts> seedShifts = {};
229  if(corrSysts) {
230  for (double systshift:{+0.5,-0.5}){
231  SystShifts tempShifts (&kMECq0ShapeSyst2017,systshift);
232  seedShifts.push_back(tempShifts);
233  }
234  }
235 
236  double minchi23 = 1E20;
237  for(double thseed:{0.45,0.55}){
238  for(int hie:{-1,1}){
239  std::cout << "\n\nFinding best fit "
240  << (thseed<0.5?"LO ":"UO ")
241  << (hie>0? "NH " : "IH ")
242  << "\n\n";
243  auto thisminchi = fit23.Fit(calc,auxShifts ,{
244  {&kFitDmSq32,{hie*fabs(calc->GetDmsq32())}},
245  {&kFitSinSqTheta23,{thseed}},
246  {&kFitDeltaInPiUnits,{0,0.5,1,1.5}}
247  }, seedShifts
248  );
249  minchi23= min(minchi23, thisminchi);
250  ResetOscCalcToDefault(calc);
251  auxShifts.ResetToNominal();
252  }//end hie
253  }//end th23
254  std::cout << "\nFound overall minchi " << minchi23 << std::endl;
255 
256  //We'll save to this file, make pretty plots later
257  TFile * outfile = new TFile(outfilename+".root","recreate");
258  //The FC also needs to know about the fitted syst values, save those in
259  //a separate file
260  TFile * outFCfile = new TFile(outFCfilename+".root","recreate");
261 
262  outfile->cd();
263  TVectorD v(1);
264  v[0] = minchi23;
265  v.Write("overall_minchi");
266 
267  //HACK HACK
268  int steps = 30;//10;// 5;
269  if(dmsqSurf){
270  std::cerr << "\n WARNING Using 20 bins for dmsq surface\n\n";
271  steps = 20;
272  }
273  //Now create all the surfaces
274  for(int hie: {-1, +1}){
275 
276  if((onlyNH && hie<0) || (onlyIH && hie>0)) continue;
277 
278  std::cout << "Starting surface " << (hie>0? "NH ": "IH") << "\n\n";
279  std::string hieStr = hie>0 ? "NH" : "IH";
280 
281  if(!th13Surf && ! dmsqSurf){
282  ResetOscCalcToDefault(calc);
283  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
284 
285  FrequentistSurface surf23(exptThis, calc,
286  &kFitDeltaInPiUnits,steps, 0, 2,
287  &kFitSinSqTheta23, steps, (nueOnly?0:0.3), (nueOnly?1:0.7),
288  {&kFitDmSq32Scaled,&kFitSinSq2Theta13},systs,
289  {},seedShifts);
290 
291  auto surf1=surf23.ToTH2(minchi23);
292  outfile->cd();
293  surf1->Write((TString)"delta_th23_"+hieStr);
294  surf23.SaveTo(outfile, (TString)"surf_delta_th23_"+hieStr);
295 
296  outFCfile->cd();
297  std::vector<TH2*> profhists = surf23.GetProfiledHists();
298  //Osc parameters are special case
299  profhists[0]->Write((hieStr+"_DmSq32").c_str());
300  profhists[1]->Write((hieStr+"_SinSq2Theta13").c_str());
301  //Also need to save a hist for each syst (new for 2017)
302  for (int i = 0; i < int(systs.size()); i++){
303  profhists[i+2]->Write((hieStr+"_"+systs[i]->ShortName()).c_str());
304  }
305  }
306  if(!th13Surf && dmsqSurf){
307  calc = DefaultOscCalc();
308  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
309 
310  FrequentistSurface surf23m(exptThis, calc,
311  &kFitSinSqTheta23, steps,0.3,0.7,
312  &kFitDmSq32Scaled, steps,
313  (hie>0?2.1:-2.8),(hie>0?2.8:-2.1),
315  {{&kFitDeltaInPiUnits,{0,0.5,1,1.5}}},seedShifts);
316  auto surf6=surf23m.ToTH2(minchi23);
317  outfile->cd();
318  surf6->Write((TString)"th23_dm32_"+(hie>0?"NH":"IH"));
319  surf23m.SaveTo(outfile, (TString)"surf_th23_dm32_"+(hie>0?"NH":"IH"));
320 
321  outFCfile->cd();
322  std::vector<TH2*> profhists = surf23m.GetProfiledHists();
323  //Osc parameters are special case
324  profhists[0]->Write((hieStr+"_SinSq2Theta13").c_str());
325  profhists[1]->Write((hieStr+"_DeltaCPpi").c_str());
326  //Also need to save a hist for each syst (new for 2017)
327  for (int i = 0; i < int(systs.size()); i++){
328  profhists[i+2]->Write((hieStr+"_"+systs[i]->ShortName()).c_str());
329  }
330 
331  }
332  if(th13Surf){
333  calc = DefaultOscCalc();
334  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
335 
336  FrequentistSurface surf13(exptThis, calc,
337  &kFitSinSq2Theta13, steps, 0, 0.35,
338  &kFitDeltaInPiUnits,steps, 0, 2,
339  {&kFitSinSqTheta23, &kFitDmSq32Scaled}, systs,
340  {{&kFitSinSqTheta23,{0.45,0.55}}}, seedShifts);
341  auto surf4 = surf13.ToTH2(minchi23);
342  outfile->cd();
343  surf4->Write((TString)"th13_delta_"+(hie>0?"NH":"IH"));
344  surf13.SaveTo(outfile, (TString)"surf_th13_delta_"+(hie>0?"NH":"IH"));
345 
346  outFCfile->cd();
347  std::vector<TH2*> profhists = surf13.GetProfiledHists();
348  //Osc parameters are special case
349  profhists[0]->Write((hieStr+"_SinSqTheta23").c_str());
350  profhists[1]->Write((hieStr+"_DmSq32").c_str());
351  //Also need to save a hist for each syst (new for 2017)
352  for (int i = 0; i < int(systs.size()); i++){
353  profhists[i+2]->Write((hieStr+"_"+systs[i]->ShortName()).c_str());
354  }
355 
356  }
357  } //end hie
358 
359  delete outfile;
360  std::cout << "\nSurfaces saved to " << outfilename << std::endl;
361  return;
362  }//end createfile
363 
364 //////////////////////////////////////////////////////////////////////
365 // Make plots
366 //////////////////////////////////////////////////////////////////////
367  else{ //plot only
368  TFile * infile = new TFile (outfilename+".root","read");
369  TFile * official =0;
370  if(corrSysts && options.Contains("joint") && fccorr && 0) {
371  official = new TFile("NOvA_nue_official_results.root","recreate");
372  TH2* axes = new TH2F("axes",";#delta_{CP}/#pi;sin^{2}#theta_{23}",100,0,2,100,0.225,0.725);
373  CenterTitles(axes);
374  axes->Write();
375  }
376 
377  auto mins =* (TVectorD*)infile->Get("overall_minchi");
378  double minchi23 = mins[0];
379 
380  std::vector <TString> surfNames;// = {"delta_th23_"};<
381  if(!th13Surf && ! dmsqSurf)surfNames = {"delta_th23_"};
382  if (dmsqSurf) surfNames.push_back("th23_dm32_");
383  if (th13Surf) surfNames.push_back("th13_delta_");
384  for(TString surfName:surfNames){
385 
386  for (int hie:{-1,+1}){
387 
388  if((onlyNH && hie<0) || (onlyIH && hie>0)) continue;
389 
390  //Setup canvas
391  if(surfName.Contains("delta_th23")) {
392  if (zoomIn) DrawContourCanvas(surfName, 0., 2., 0.225, 0.725);
393  else DrawContourCanvas(surfName, 0., 2., 0., 1.);
394  }
395  if(surfName.Contains("th23_dm32")) {
396  if (hie>0) DrawContourCanvas(surfName, 0.3, 0.7, 2.05, 2.75);
397  else DrawContourCanvas(surfName, 0.3, 0.7, -2.95, -2.25);
398  }
399  if(surfName.Contains("th13_delta")) {
400  if (zoomIn) DrawContourCanvas(surfName, 0., 0.38, 0, 2);
401  else DrawContourCanvas(surfName, 0., 0.5, 0, 2);
402  }
403  //Load the surface; fill and FC corrections if available
404  auto surf = *ISurface::LoadFrom(infile, "surf_" + surfName + (hie > 0 ? "NH" : "IH"));
405 
406  TH2 * surf_1Sigma = Gaussian68Percent2D(surf);
407  TH2 * surf_2Sigma = Gaussian2Sigma2D(surf);
408  TH2 * surf_3Sigma = Gaussian3Sigma2D(surf);
409 
410  if(fccorr){
411  if(!joint || !corrSysts || decomp != "combo" || dmsqSurf || th13Surf)
412  { std::cerr << "\nWe don't have FC\n\n"; exit(1);}
413 
414  TString fcFolder="/nova/ana/nu_e_ana/Ana2017/FC/";
415  TString fcFileName ((TString)"FCCol_"+(hie>0?"NH":"IH")+"_surf_ssth23delta.root");
416  FCSurface *fcXH = 0;
417  auto fccol = FCCollection::LoadFromFile((fcFolder + fcFileName).Data()).release();
418  fcXH = new FCSurface(30,0,2,30,0.3,0.7);
419  fcXH->Add(*fccol);
420  surf_1Sigma = fcXH->SurfaceForSignificance(0.6827);
421  surf_2Sigma = fcXH->SurfaceForSignificance(0.9545);
422  surf_3Sigma = fcXH->SurfaceForSignificance(0.9973);
423  if(smooth){
424  surf_1Sigma->Smooth();
425  surf_2Sigma->Smooth();
426  surf_3Sigma->Smooth();
427  }
428  }
429 
430  Color_t k2SigmaConfidenceColor = ((hie>0)? k2SigmaConfidenceColorNH:k2SigmaConfidenceColorIH);
431  Int_t kFillColor1 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.70);
432  Int_t kFillColor2 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.35);
433  Int_t kFillColor3 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.20);
434  Int_t kDarkColor = kGray+3;
435 
436  if(fillContour){
437  auto g1 = surf.GetGraphs(surf_1Sigma,minchi23);
438  auto g2 = surf.GetGraphs(surf_2Sigma,minchi23);
439  auto g3 = surf.GetGraphs(surf_3Sigma,minchi23);
440  FillGraphs(g1, g2, g3,
441  kFillColor1, kFillColor2, kFillColor3,
442  kDarkColor, suffix + surfName + (hie>0?"NH":"IH") + (fccorr?"_fccorr":""));
443  }
444  else{
445  surf.DrawContour(surf_1Sigma, kSolid, kFillColor1,minchi23);
446  surf.DrawContour(surf_2Sigma, k2SigmaConfidenceStyle, kFillColor2,minchi23);
447  surf.DrawContour(surf_3Sigma, kSolid, kFillColor3,minchi23);
448  }
449 
450  gPad->RedrawAxis();
451 
452  bool bestfit = 0; //todo
453  auto leg = ContourLegend(hie, fillContour, fccorr,
454  kFillColor1,kFillColor2,kFillColor3,kDarkColor,
455  bestfit);
456  if(th13Surf) {
457  leg->SetFillStyle(1001);
458  leg->SetNColumns(1);
459  leg->SetX1(0.76); leg->SetX2(0.89);
460  leg->SetY1(0.52); leg->SetY2(0.85);
461  leg->SetHeader("No FC");
462  }
463  leg->Draw();
464  DrawHieTag(hie,false);
465 
466  for(TString ext: {".pdf",".root"}){
467  gPad->Print(outdir + "contour_" + suffix + "_" +
468  surfName +
469  (hie>0?"NH":"IH") + (fccorr?"_fccorr":"") +
470  ((fccorr && smooth)?"_smooth":"") +
471  (zoomIn?"_zoom":"") +
472  ext);
473  }//end ext
474 
475  if(debug && !fccorr){
476  auto hists = surf.GetProfiledHists();
477  int nSysts = corrSysts ? 30 : 0;
478  for (int i = 0; i < nSysts; ++i){
479  hists.push_back((TH2*) infile->Get((TString) "surf_" + surfName +
480  (hie > 0 ? "NH" : "IH") +
481  "/profHists/hist" + std::to_string(i)));
482  }
483  hists.push_back(surf.ToTH2(minchi23));
484 
485  for(uint i = 0; i < hists.size(); ++i){
486  if(!hists[i]) continue;
487  if(TString(hists[i]->GetTitle()).Contains("delta")) {
489  hists[i]->GetZaxis()->SetRangeUser(0,2);
490  }
491  else SetPaletteBlueRedWhite();
492  if( i > 3) {
493  double minz = hists[i]->GetBinContent(hists[i]->GetMinimumBin());
494  double maxz = hists[i]->GetBinContent(hists[i]->GetMaximumBin());
495  double lim = std::max(fabs(minz), fabs(maxz));
496  hists[i]->GetZaxis()->SetRangeUser(-lim,+lim);
497  }
498  if(i==(hists.size()-1)) {
500  hists[i]->GetZaxis()->SetRangeUser(0.1,20);
501  gPad->SetLogz();
502  }
503 
504  hists[i]->Draw("colz");
505  hists[i]->GetYaxis()->SetTitleOffset(0.8);
506  gPad->SetRightMargin(0.14);
507  surf.DrawContour(surf_1Sigma, kSolid, kDarkColor,minchi23);
508  surf.DrawContour(surf_2Sigma, kSolid, kDarkColor,minchi23);
509  surf.DrawContour(surf_3Sigma, kSolid, kDarkColor,minchi23);
510  gPad->Print(outdir + "debug_contour_" + suffix + "_" +
511  surfName +
512  (hie > 0 ? "NH":"IH") + (fccorr?"_fccorr":"") +
513  (zoomIn ? "_zoom" : "") +
514  (std::to_string(i)) + ".pdf");
515  }//end profiles
516  }//end debug
517  }//end hie
518  }//end surf
519  } // end plot only
520 }//end all
521 
522 
523 void DrawHieTag(int hie,bool th13, bool vert){
524  if(hie>0){
525  TLatex * ltxNH = new TLatex(0.16,0.89,"#bf{NH}");
526  ltxNH->SetTextColor(kPrimColorNH);
527  ltxNH->SetNDC();
528  ltxNH->SetTextAlign(22);
529  ltxNH->SetTextSize(kBlessedLabelFontSize);
530  if(th13) ltxNH->Draw();
531  else if(!vert) ltxNH->DrawLatex(.85,.20,"#bf{NH}");
532  else ltxNH->DrawLatex(.85,.85,"#bf{NH}");
533  }
534  else{
535  TLatex * ltxIH = new TLatex (0.16,0.46,"#bf{IH}");
536  ltxIH->SetTextColor(kPrimColorIH);
537  ltxIH->SetNDC();
538  ltxIH->SetTextAlign(22);
539  ltxIH->SetTextSize(kBlessedLabelFontSize);
540  if(th13) ltxIH->Draw();
541  else if(!vert) ltxIH->DrawLatex(.85,.20,"#bf{IH}");
542  else ltxIH->DrawLatex(.85,.85,"#bf{IH}");
543  }
544 }
545 
546 TLegend * ContourLegend(int hie, bool fillContour, bool fccorr,
547  Int_t kFillColor1,Int_t kFillColor2, Int_t kFillColor3,Int_t kDarkColor,
548  bool bestFit){
549  TLegend * leg = new TLegend(0.16,fccorr?0.19:0.17,0.65,fccorr?0.26:0.28);
550  leg->SetTextSize(kBlessedLabelFontSize);
551  leg->SetFillStyle(0);
552  leg->SetMargin(1.3*leg->GetMargin());
553  leg->SetNColumns(3);
554 
555  TGraph * dummy = new TGraph;
556  dummy->SetLineWidth(3);
557 
558  if(fillContour){
559  dummy->SetLineColor(kDarkColor);
560  dummy->SetFillColor(kFillColor1);
561  leg->AddEntry(dummy->Clone(),"1 #sigma","f");
562  dummy->SetFillColor(kFillColor2);
563  leg->AddEntry(dummy->Clone(),"2 #sigma","f");
564  dummy->SetFillColor(kFillColor3);
565  leg->AddEntry(dummy->Clone(),"3 #sigma","f");
566  }
567  else{
568  leg->SetMargin(1.4*leg->GetMargin());
569  dummy->SetLineColor(kFillColor1);
570  dummy->SetLineStyle(kSolid);
571  leg->AddEntry(dummy->Clone(),"1 #sigma","l");
572  dummy->SetLineStyle(k2SigmaConfidenceStyle);
573  dummy->SetLineColor(kFillColor2);
574  leg->AddEntry(dummy->Clone(),"2 #sigma","l");
575  dummy->SetLineStyle(kSolid);
576  dummy->SetLineColor(kFillColor3);
577  leg->AddEntry(dummy->Clone(),"3 #sigma","l");
578  }
579  if(bestFit){
580  leg->SetNColumns(4);
581 // leg->SetColumnSeparation(0.08);
582  leg->SetX2(0.75);
583  dummy->SetMarkerStyle(kFullCircle);
584  dummy->SetMarkerColor(kDarkColor);
585  dummy->SetLineColor(kDarkColor);
586  dummy->SetLineStyle(7);
587  if(hie>0) leg->AddEntry(dummy->Clone(),"Best Fit","p");
588  if(hie<0) leg->AddEntry((TObject*)0, "#color[0]{Best Fit}", "");
589  dummy->SetLineStyle(1);
590  }
591  if(!fccorr) leg->SetHeader("No Feldman-Cousins");
592  return leg;
593 }
594 
595 void FillGraphs(std::vector<TGraph*> g1,std::vector<TGraph*> g2,std::vector<TGraph*> g3,
596  const Int_t kFillColor1, const Int_t kFillColor2, const Int_t kFillColor3,
597  const Int_t kDarkColor, const TString surftype)
598 {
599  bool isJoint = surftype.Contains("joint");
600  bool isNH = surftype.Contains("NH");
601  bool isSysts = surftype.Contains("systs");
602  bool isFccorr = surftype.Contains("fccorr");
603 
604  if (surftype.Contains("delta_th23")){
605  if(isJoint && isNH ){
606  JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f");
607  JoinGraphs(g2[0], g2[1], kWhite)->Draw("f");
608  JoinGraphs(g2[0], g2[1], kFillColor2)->Draw("f");
609  if(isSysts && !isFccorr) JoinGraphs(g1[0], g1[1], kWhite)->Draw("f");
610  if(isSysts && !isFccorr) JoinGraphs(g1[0], g1[1], kFillColor1)->Draw("f");
611  if(!isSysts || isFccorr) for (auto & g:g1) { g->SetFillColor(kWhite); g->DrawClone("f");}
612  if(!isSysts || isFccorr) for (auto & g:g1) { g->SetFillColor(kFillColor1); g->Draw("f");}
613  }
614 
615  if(isJoint && !isNH){
616  for (auto &g:g3) {
617  g->SetFillColor(kFillColor3); g->Draw("f");}
618  for (auto &g:g2) {
619  g->SetFillColor(kWhite); g->DrawClone("f");
620  g->SetFillColor(kFillColor2); g->Draw("f");
621  }
622  }
623  }
624 
625  if (surftype.Contains("th23_dm32")){
626  for (auto &g:g3) {
627  g->SetFillColor(kFillColor3); g->Draw("f");}
628  for (auto &g:g2) {
629  g->SetFillColor(kWhite); g->DrawClone("f");
630  g->SetFillColor(kFillColor2); g->Draw("f");
631  }
632  for (auto &g:g1) {
633  g->SetFillColor(kWhite); g->DrawClone("f");
634  g->SetFillColor(kFillColor1); g->Draw("f");
635  }
636  }
637 
638  if (surftype.Contains("th13_delta")) {
639  JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f");
640  JoinGraphs(g2[0], g2[1], kWhite)->Draw("f");
641  JoinGraphs(g2[0], g2[1], kFillColor2)->Draw("f");
642  JoinGraphs(g1[0], g1[1], kWhite)->Draw("f");
643  JoinGraphs(g1[0], g1[1], kFillColor1)->Draw("f");
644  }
645 
646 
647  for (auto & gs:{g3,g2,g1}){
648  for (auto & g:gs) {
649  g->SetLineColor(kDarkColor);
650  g->SetLineWidth(3);
651  g->Draw("l");
652  }
653  }
654 
655 }
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
std::vector< std::pair< TH1D *, double > > GetNumuCosmics2017(const int nq=4, bool GetFromUPS=false)
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
const Color_t kCosmicBackgroundColor
Definition: Style.h:41
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:18
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
double maxy
const Dmsq32Constraint kDmsq32ConstraintPDG2017(2.45e-3, 0.05e-3, 2.52e-3, 0.05e-3)
void XAxisDeltaCPLabels(TH1 *axes, bool t2kunits)
Label the x-axis with fractions of pi.
TH2 * Gaussian68Percent2D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 2D in gaussian approximation.
OStream cerr
Definition: OStream.cxx:7
const IPrediction * GetNuePrediction2017(std::string decomp, osc::IOscCalc *calc, bool corrSysts, bool GetFromUPS=false)
const double kAna2017POT
Definition: Exposures.h:174
const Color_t k2SigmaConfidenceColorNH
Definition: Style.h:77
TH2 * Gaussian2Sigma2D(const FrequentistSurface &s)
Up-value surface for 2 sigma confidence in 2D in gaussian approximation.
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1481
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 XML_Char const XML_Char * data
Definition: expat.h:268
Log-likelihood scan across two parameters.
string outfilename
knobs that need extra care
std::pair< TH1D *, double > GetNueCosmics2017(bool GetFromUPS=false)
void SetPaletteBlueRedWhite()
Definition: Style.cxx:7
void DrawHieTag(int hie, bool th13=true, bool vert=false)
std::vector< const IPrediction * > GetNumuPredictions2017(const int nq=4, bool useSysts=true, bool GetFromUPS=false)
string infile
std::vector< Spectrum * > GetNumuData2017(const int nq=4)
void SetPaletteWhiteBlueDark()
Definition: Style.cxx:34
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
Spectrum * GetNueData2017()
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.
void SetPaletteBlueRedCyclic()
Definition: Style.cxx:19
virtual T GetDmsq32() const
Definition: IOscCalc.h:86
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 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
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
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.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const Style_t k2SigmaConfidenceStyle
Definition: Style.h:71
const NOvARwgtSyst kMECq0ShapeSyst2017("MECq0Shape","MEC q_{0} shape", novarwgt::kMECq0ShapeSyst2017)
Definition: MECSysts.h:41
void ResetToNominal()
Definition: SystShifts.cxx:141
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 ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
exit(0)
assert(nhit_max >=nhit_nbins)
c1
Definition: demo5.py:24
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
TH1 * DataMCComparison(const Spectrum &data, const Spectrum &mc, EBinType bintype)
Definition: Plots.cxx:35
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
void joint_fit_2017_contours(std::string decomp="prop", bool createFile=false, bool corrSysts=false, TString options="joint_fakeNHLO", bool dmsqSurf=false, bool th13Surf=false, bool fccorr=false, bool zoomIn=true, bool fillContour=false, bool smooth=true)
TH2 * SurfaceForSignificance(double sig)
Definition: FCSurface.cxx:103
Float_t e
Definition: plot.C:35
const Float_t kBlessedTitleFontSize
Definition: Style.h:89
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)
double th13
Definition: runWimpSim.h:98
static std::unique_ptr< FCCollection > LoadFromFile(const std::string &wildcard)
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:30
const std::string outdir
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
unsigned int uint