plot_joint_fit_2020_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 
20 #include "../joint_fit_2020_loader_tools.h"
21 #include "../joint_fit_2020_datarelease_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 plot_joint_fit_2020_contours(bool corrSysts = false,
67  bool runOnGrid = false,
68  TString options="joint_fake2019_both_onlyIH",
69  bool dmsqSurf = false,
70  bool th13Surf = false,
71  bool fccorr=false,
72  bool zoomIn = true,
73  bool fillContour = false,
74  bool smooth = true,
75  bool fc_overlay = false)
76 {
77  bool debug = 1;
78  bool nueOnly = options.Contains("nueOnly");
79  bool numuOnly = options.Contains("numuOnly");
80  bool joint = options.Contains("joint");
81  assert (nueOnly || numuOnly || joint);
82 
83  bool FHCOnly = options.Contains("FHCOnly");
84  bool RHCOnly = options.Contains("RHCOnly");
85  bool both = options.Contains("both");
86  assert (FHCOnly || RHCOnly || both);
87 
88  bool fake2019 = options.Contains("fake2019");
89  bool realData = options.Contains("realData");
90 
91  bool onlyNH = options.Contains("onlyNH");
92  bool onlyIH = options.Contains("onlyIH");
93 
94  auto suffix = options;
95  if(dmsqSurf) suffix+= "_dmsq";
96  if(th13Surf) suffix+= "_th13";
97  if(corrSysts) suffix+="_systs";
98  else suffix+="_stats";
99 
100  TString outdir = "/nova/ana/3flavor/Ana2020/Fits/RealData/WithSeeds/";
101  if (runOnGrid) outdir = "./";
102  TString outFCfilename (outdir +"contours_FCInput_2020_" + suffix);
103  TString outfilename (outdir + "hist_contours_2020_" + suffix);
104 
105 
106  //////////////////////////////////////////////////////////////////////
107  // Make plots
108  //////////////////////////////////////////////////////////////////////
109 
110  /*TString infilename = "/nova/ana/nu_e_ana/Ana2019/Results/RHC_and_FHC/contours/";
111  if(th13Surf) infilename += "th13_delta/";
112  else if (dmsqSurf) infilename += "th23_dmsq/";
113  else infilename += "delta_th23/";
114  infilename += corrSysts?"syst/":"stat/";*/
115  TString infilename = outdir;
116  infilename += "hist_contours_2020_";
117  infilename += options;
118  if (dmsqSurf) infilename += "_dmsq";
119  infilename += corrSysts?"_systs.root":"_stats.root";
120 
121  TFile * infile = new TFile (infilename,"read");
122  TFile * official =0;
123  TString officialpreffix = "./";
124  bool makedatarelease = true;
125 
126  if(!makedatarelease && corrSysts && options.Contains("joint") && fccorr && 0) {
127  official = new TFile("NOvA_nue_official_results.root","recreate");
128  TH2* axes = new TH2F("axes",";#delta_{CP}/#pi;sin^{2}#theta_{23}",100,0,2,100,0.225,0.725);
129  //CenterTitles(axes);
130  axes->Write();
131  }
132  else{
133  if (makedatarelease){
134  official = MakeDataReleaseFileContour(officialpreffix, dmsqSurf);
135  }
136  }
137 
138  auto mins =* (TVectorD*)infile->Get("overall_minchi");
139  double minchi23 = mins[0];
140 
141  std::vector <TString> surfNames;// = {"delta_th23_"};<
142  if(!th13Surf && ! dmsqSurf)surfNames = {"delta_th23_"};
143  if (dmsqSurf) surfNames.push_back("th23_dm32_");
144  if (th13Surf) surfNames.push_back("th13_delta_");
145 
146  for(TString surfName:surfNames){
147  for (int hie:{-1,+1}){
148 
149  if((onlyNH && hie<0) || (onlyIH && hie>0)) continue;
150 
151  //Setup canvas
152  if(surfName.Contains("delta_th23")) {
153  if (zoomIn) DrawContourCanvas(surfName, 0., 2., (RHCOnly?0:0.275), (RHCOnly?1.0:0.725));
154  else DrawContourCanvas(surfName, 0., 2., 0., 1.);
155  }
156  if(surfName.Contains("th23_dm32")) {
157  if (hie>0) DrawContourCanvas(surfName, ((RHCOnly)?0.2:0.3), ((RHCOnly)?0.8:0.7), (RHCOnly?2.0:2.075), (RHCOnly?3.3:2.725));
158  else DrawContourCanvas(surfName, ((RHCOnly)?0.2:0.3), ((RHCOnly)?0.8:0.7), (RHCOnly?-3.0:-2.8), (RHCOnly?-2.0:-2.2));
159  }
160  if(surfName.Contains("th13_delta")) {
161  if (zoomIn) DrawContourCanvas(surfName, 0., 0.38, 0, 2);
162  else DrawContourCanvas(surfName, 0., 1.0, 0, 2);
163  }
164  //Load the surface; fill and FC corrections if available
165  auto surf = *FrequentistSurface::LoadFrom(infile, "surf_"+surfName+(hie>0?"NH":"IH"));
166  TH2 * hsurf = (TH2*)infile->Get(surfName+(hie>0?"NH":"IH"));
167  TH2 * surf_1Sigma = Gaussian68Percent2D(surf);
168  TH2 * surf_2Sigma = Gaussian2Sigma2D(surf);
169  TH2 * surf_3Sigma = Gaussian3Sigma2D(surf);
170  TH2 * fc_surf_1Sigma = new TH2F();
171  TH2 * fc_surf_2Sigma = new TH2F();
172  TH2 * fc_surf_3Sigma = new TH2F();
173  TH2 * fc_surf_90CL = new TH2F();
174  if(fccorr){
175  // if(!joint || !corrSysts || decomp != "combo" || dmsqSurf || th13Surf)
176  // { std::cerr << "\nWe don't have FC\n\n"; exit(1);}
177 
178  TString fcFolder="/nova/ana/nu_e_ana/Ana2020/FC/";
179  //TString fcFolder = outdir;
180  TString fcFileName;
181  if(dmsqSurf) fcFileName = "contours/ssth23dmsq/ssth23dmsq_";
182  else if(!dmsqSurf && !th13Surf)
183  fcFileName = "contours/deltassth23/deltassth23_";
184  fcFileName += hie > 0? "nh":"ih";
185  fcFileName += ".root";
186 
187  FCSurface *fcXH = 0;
188  auto fccol = FCCollection::LoadFromFile((fcFolder + fcFileName).Data()).release();
189  int NX = hsurf->GetXaxis()->GetNbins();
190  int NY = hsurf->GetYaxis()->GetNbins();
191  double minX = hsurf->GetXaxis()->GetBinLowEdge(1);
192  double maxX = hsurf->GetXaxis()->GetBinUpEdge(NX);
193  double minY = hsurf->GetYaxis()->GetBinLowEdge(1);
194  double maxY = hsurf->GetYaxis()->GetBinUpEdge(NY);
195 
196  std::cout << "Making FCSurface: \n";
197  std::cout << "-------------------------------\n";
198  std::cout << "NX: " << NX;
199  std::cout << "\tNY: " << NY;
200  std::cout << "\tminX: " << minX;
201  std::cout << "\tmaxX: " << maxX;
202  std::cout << "\tminY: " << minY;
203  std::cout << "\tmaxY: " << maxY;
204  std::cout << "\n------------------------------\n";
205 
206  fcXH = new FCSurface(NX, minX, maxX, NY, minY, maxY);
207  std::string plot_name = "deltassth23";
208  if(dmsqSurf) plot_name = "ssth23dmsq32";
209  if(th13Surf) plot_name = "deltath13";
210 
211  fcXH->Add(*fccol, plot_name);
212  fc_surf_1Sigma = fcXH->SurfaceForSignificance(0.6827);
213  fc_surf_2Sigma = fcXH->SurfaceForSignificance(0.9545);
214  fc_surf_3Sigma = fcXH->SurfaceForSignificance(0.9973);
215  fc_surf_90CL = fcXH->SurfaceForSignificance(0.90);
216  if(smooth){
217  fc_surf_1Sigma->Smooth();
218  fc_surf_2Sigma->Smooth();
219  fc_surf_90CL->Smooth();
220  if(dmsqSurf || th13Surf)
221  fc_surf_3Sigma->Smooth();
222  else
223  CylindricalSmoothing(fc_surf_3Sigma);
224  }
225  }
226  Color_t k2SigmaConfidenceColor = ((hie>0)? cnh:cih);
227  Int_t kFillColor1 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.99);
228  Int_t kFillColor2 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.5);
229  Int_t kFillColor3 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.2);
230  Int_t kFCFillColor1 = TColor::GetColorTransparent(k2SigmaConfidenceColor, fc_overlay? 0.65: 0.75);
231  Int_t kFCFillColor2 = TColor::GetColorTransparent(k2SigmaConfidenceColor, fc_overlay? 0.4: 0.38);
232  Int_t kFCFillColor3 = TColor::GetColorTransparent(k2SigmaConfidenceColor, fc_overlay? 0.18: 0.18);
233  Int_t kDarkColor = kGray+3;
234 
235  if(fillContour){
236  if(!fccorr){
237  auto g1 = surf.GetGraphs(surf_1Sigma,minchi23);
238  auto g2 = surf.GetGraphs(surf_2Sigma,minchi23);
239  auto g3 = surf.GetGraphs(surf_3Sigma,minchi23);
240  FillGraphs(g1, g2, g3,
241  kFillColor1, kFillColor2, kFillColor3,
242  kDarkColor, suffix + surfName + (hie>0?"NH":"IH"));
243  }
244  else {
245  auto fc_g1 = surf.GetGraphs(fc_surf_1Sigma,minchi23);
246  auto fc_g2 = surf.GetGraphs(fc_surf_2Sigma,minchi23);
247  auto fc_g3 = surf.GetGraphs(fc_surf_3Sigma,minchi23);
248  FillGraphs(fc_g1, fc_g2, fc_g3,
249  kFillColor1, kFillColor2, kFillColor3,
250  kDarkColor, suffix + surfName + (hie>0?"NH":"IH")+ "fccorr");
251  }
252  }
253  else{
254  cout<<"minchi from file is "<<minchi23<<endl;
255  if(!fccorr || fc_overlay){
256 
257  if(fc_overlay) {
258  auto g1 = surf.GetGraphs(surf_1Sigma,minchi23);
259  auto g2 = surf.GetGraphs(surf_2Sigma,minchi23);
260  auto g3 = surf.GetGraphs(surf_3Sigma,minchi23);
261  FillGraphs(g1, g2, g3,
262  kFillColor1, kFillColor2, kFillColor3,
263  kDarkColor, suffix + surfName + (hie>0?"NH":"IH"));
264  }
265  surf.DrawContour(surf_1Sigma, kSolid, kFillColor1,minchi23);
266  surf.DrawContour(surf_2Sigma, k2SigmaConfidenceStyle, kFillColor2,minchi23);
267  surf.DrawContour(surf_3Sigma, kSolid, kFillColor3,minchi23);
268  }
269  if(fccorr) {
270  surf.DrawContour(fc_surf_1Sigma, kSolid, kFCFillColor1, minchi23);
271  surf.DrawContour(fc_surf_2Sigma, k2SigmaConfidenceStyle, kFCFillColor2, minchi23);
272  surf.DrawContour(fc_surf_3Sigma, kSolid, kFCFillColor3, minchi23);
273  }
274  }
275 
276  gPad->RedrawAxis();
277 
278  bool bestfit = 1;
279  double delta, th23, dmsq32, th13;
280  TMarker bf;
281  if(bestfit && both && corrSysts && hie>0){
282  TFile * bffile = new TFile("/nova/ana/3flavor/Ana2020/Fits/RealData/WithSeeds/bestfits_joint_realData_both_systs.root","READ");
283  auto calc_tmp = LoadFrom<osc::IOscCalc>(bffile, "NH_UO_calc").release();
284  auto calc = (osc::IOscCalcAdjustable*)calc_tmp->Copy();
285  if(!th13Surf){
288  dmsq32 = kFitDmSq32Scaled.GetValue(calc);
290  }
291  else{delta = 0.678516; th13 = 0.0957894; cout<<"NB! Old 2019 values"<<endl;} /// NB! Old 2019 values
292  bf.SetMarkerColor(kDarkColor);
293  // bf.SetMarkerStyle(kFullStar);
294  bf.SetMarkerStyle(43);
295  bf.SetMarkerSize(2);
296  if(!dmsqSurf && ! th13Surf){
297  bf.SetX(delta);
298  bf.SetY(th23);
299  bf.DrawClone();
300  }
301  if(dmsqSurf){
302  bf.SetX(th23);
303  bf.SetY(dmsq32);
304  bf.DrawClone();
305  }
306  if(th13Surf){
307  bf.SetX(th13);
308  bf.SetY(delta);
309  bf.DrawClone();
310  }
311  }
312  auto leg = ContourLegend(hie, fillContour, fccorr,
313  kFillColor1,kFillColor2,kFillColor3,kDarkColor,
314  bestfit && both);
315 
316  if(th13Surf) {
317  DrawReactorConstraint(0.082,0.004);
318  leg->SetFillStyle(1001);
319  leg->SetNColumns(1);
320  leg->SetX1(0.7); leg->SetX2(0.8);
321  leg->SetY1(0.4); leg->SetY2(0.85);
322  if(!fccorr) leg->SetHeader("No FC");
323 
324  TGraph * dummy = new TGraph;
325  dummy->SetLineWidth(2);
326  dummy->SetLineColor(kGray+2);
327  dummy->SetFillColor(kGray);
328  leg->AddEntry(dummy,"#splitline{Reactor}{68%C.L.}","f");
329  }
330  leg->Draw();
331  Preliminary();
332  DrawHieTag(hie,false);
333 
334 
335  for(TString ext: {".png", ".pdf",".root"}){
336  gPad->Print("./contour_" + suffix + "_" +
337  surfName +
338  (hie>0?"NH":"IH") + (fccorr?"_fccorr":"") +
339  ((fccorr && smooth)?"_smooth":"") +
340  (zoomIn?"_zoom":"") +
341  (fc_overlay?"_overlay":"") +
342  ext);
343  }//end ext
344 
345  if(official) {
346  auto fc_g1 = surf.GetGraphs(fc_surf_1Sigma,minchi23);
347  auto fc_g2 = surf.GetGraphs(fc_surf_2Sigma,minchi23);
348  auto fc_g3 = surf.GetGraphs(fc_surf_3Sigma,minchi23);
349  auto fc_g9 = surf.GetGraphs(fc_surf_90CL,minchi23);
350  GraphsToDataRelease(official, fc_g1,"1sigma", hie>0?"NH":"IH");
351  GraphsToDataRelease(official, fc_g2,"2sigma", hie>0?"NH":"IH");
352  GraphsToDataRelease(official, fc_g3,"3sigma", hie>0?"NH":"IH");
353  GraphsToDataRelease(official, fc_g9,"90CL", hie>0?"NH":"IH");
354  if(hie>0) //gah
355  BestFitToDataRelease(official, bf, "best_fit",hie>0?"NH":"IH");
356  }
357 
358  debug = false;
359  if(debug && !fccorr){
360  auto hists = surf.GetProfiledHists();
361  int nSysts = corrSysts ? 70 : 0;
362  for (int i = 0; i < nSysts; ++i){
363  hists.push_back((TH2*) infile->Get((TString) "surf_" + surfName +
364  (hie > 0 ? "NH" : "IH") +
365  "/profHists/hist" + std::to_string(i)));
366  }
367  hists.push_back(surf.ToTH2(minchi23));
368 
369  for(uint i = 0; i < hists.size(); ++i){
370  if(!hists[i]) continue;
371  if(TString(hists[i]->GetTitle()).Contains("delta")) {
373  hists[i]->GetZaxis()->SetRangeUser(0,2);
374  }
375  else gStyle->SetPalette(87);
376  if( i > 3) {
377  double minz = hists[i]->GetBinContent(hists[i]->GetMinimumBin());
378  double maxz = hists[i]->GetBinContent(hists[i]->GetMaximumBin());
379  double lim = std::max(fabs(minz), fabs(maxz));
380  hists[i]->GetZaxis()->SetRangeUser(-lim,+lim);
381  }
382  if(i==(hists.size()-1)) {
383  gStyle->SetPalette(56);
384  hists[i]->GetZaxis()->SetRangeUser(0.,11.83);
385  }
386 
387  hists[i]->Draw("colz");
388  hists[i]->GetYaxis()->SetTitleOffset(0.8);
389  gPad->SetRightMargin(0.14);
390  surf.DrawContour(surf_1Sigma, kSolid, kDarkColor,minchi23);
391  surf.DrawContour(surf_2Sigma, kSolid, kDarkColor,minchi23);
392  surf.DrawContour(surf_3Sigma, kSolid, kDarkColor,minchi23);
393  gPad->Print("./debug_contour_" + suffix + "_" +
394  surfName +
395  (hie > 0 ? "NH":"IH") + (fccorr?"_fccorr":"") +
396  (zoomIn ? "_zoom" : "") +
397  (std::to_string(i)) + ".pdf");
398 
399  }//end profiles
400  }//end debug
401  }//end hie
402  }//end surf
403 }//end all
404 
405 
406 void DrawHieTag(int hie,bool th13, bool vert){
407  if(hie>0){
408  TLatex * ltxNH = new TLatex(0.16,0.89,"#bf{NH}");
409  ltxNH->SetTextColor(cnh);
410  ltxNH->SetNDC();
411  ltxNH->SetTextAlign(22);
412  ltxNH->SetTextSize(kBlessedLabelFontSize);
413  if(th13) ltxNH->Draw();
414  else if(!vert) ltxNH->DrawLatex(.85,.20,"#bf{NH}");
415  else ltxNH->DrawLatex(.85,.85,"#bf{NH}");
416  }
417  else{
418  TLatex * ltxIH = new TLatex (0.16,0.46,"#bf{IH}");
419  ltxIH->SetTextColor(cih_dark);
420  ltxIH->SetNDC();
421  ltxIH->SetTextAlign(22);
422  ltxIH->SetTextSize(kBlessedLabelFontSize);
423  if(th13) ltxIH->Draw();
424  else if(!vert) ltxIH->DrawLatex(.85,.20,"#bf{IH}");
425  else ltxIH->DrawLatex(.85,.85,"#bf{IH}");
426  }
427 }
428 
429 
430 void FillGraphs(std::vector<TGraph*> g1,std::vector<TGraph*> g2,std::vector<TGraph*> g3,
431  const Int_t kFillColor1, const Int_t kFillColor2, const Int_t kFillColor3,
432  const Int_t kDarkColor, const TString surftype)
433 {
434  bool isJoint = surftype.Contains("joint");
435  bool isBoth = surftype.Contains("both");
436  bool isRHC = surftype.Contains("RHCOnly");
437  bool isNH = surftype.Contains("NH");
438  bool isSysts = surftype.Contains("systs");
439  bool isFccorr = surftype.Contains("fccorr");
440 
441  if (surftype.Contains("delta_th23")){
442  cout<<"delta"<<endl;
443  if(isJoint && isNH ){
444  cout<<"joint nh"<<endl;
445  if(isRHC){
446  JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f");
447  JoinGraphs(g2[0], g2[1], kWhite)->Draw("f");
448  JoinGraphs(g2[1], g2[2], kWhite)->Draw("f");
449  JoinGraphs(g2[0], g2[1], kFillColor2)->Draw("f same");
450  JoinGraphs(g2[1], g2[2], kFillColor2)->Draw("f same");
451  }
452  else{
453  JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f");
454  JoinGraphs(g3[0], g3[1], kWhite)->Draw("l");
455  JoinGraphs(g2[0], g2[1], kWhite)->Draw("f");
456  JoinGraphs(g2[0], g2[1], kFillColor2)->Draw("f");
457  JoinGraphs(g2[0], g2[1], kWhite)->Draw("l");
458  //JoinGraphs(g1[0], g1[1], kFillColor1)->Draw("c f");
459  //JoinGraphs(g1[0], g1[1], kWhite)->Draw("c");
460  if(isBoth){
461  for (auto & g:g1) { g->SetFillColor(kWhite); g->DrawClone("c f");}
462  for (auto & g:g1) { g->SetFillColor(kFillColor1); g->Draw("c f");}
463  for (auto & g:g1) { g->Draw("c");}
464  }
465  else{
466  JoinGraphs(g1[0], g1[1], kFillColor1)->Draw("c f");
467  JoinGraphs(g1[0], g1[1], kWhite)->Draw("c");
468  }
469  }
470  }
471 
472  if(isJoint && !isNH){
473  if(isRHC) JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f");
474  else{
475  //JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f c");
476  g3[1]->SetFillColor(kFillColor3); g3[1]->Draw("f c");
477  g3[0]->SetFillColor(kFillColor3); g3[0]->Draw("f c");
478  for (auto &g:g3) {
479  g->SetLineWidth(3);
480  g->Draw("c");
481  }
482  }
483  for (auto &g:g2) {
484  g->SetFillColor(kWhite); g->DrawClone("f c");
485  g->SetFillColor(kFillColor2); g->Draw("f c");
486  g->SetLineColor(kDarkColor);
487  g->SetLineWidth(3);
488  g->Draw("c");
489  }
490  for (auto &g:g1) {
491  g->SetFillColor(kWhite); g->DrawClone("f c");
492  g->SetFillColor(kFillColor1); g->Draw("f c");
493  g->SetLineColor(kDarkColor);
494  g->SetLineWidth(3);
495  g->Draw("c");
496  }
497  for(auto & gs:{g3, g2, g1}){
498  for(auto & g:gs) {
499  g->SetLineColor(kDarkColor);
500  g->SetLineWidth(3);
501  g->Draw("c");
502  }
503  }
504  }
505  }
506 
507  if (surftype.Contains("th23_dm32")){
508  for (auto &g:g3) {
509  g->SetFillColor(kFillColor3); g->Draw("cf");
510  }
511  for (auto &g:g2) {
512  g->SetFillColor(kWhite); g->DrawClone("cf");
513  g->SetFillColor(kFillColor2); g->Draw("cf");
514  }
515  for (auto &g:g1) {
516  g->SetFillColor(kWhite); g->DrawClone("cf");
517  g->SetFillColor(kFillColor1); g->Draw("cf");
518  }
519  for (auto & gs:{g3,g2,g1}){
520  for (auto & g:gs) {
521  g->SetLineColor(kDarkColor);
522  g->SetLineWidth(3);
523  g->Draw("c");
524  }
525  }
526 
527  }
528 
529  if (surftype.Contains("th13_delta")) {
530  if(isBoth) {
531  JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f");
532  if(isNH) {
533  JoinGraphs(g2[0], g2[1], kWhite)->Draw("f");
534  JoinGraphs(g2[0], g2[1], kFillColor2)->Draw("f");
535  }
536  else{
537  for (auto & g:g2) { g->SetFillColor(kWhite); g->DrawClone("f");}
538  for (auto & g:g2) { g->SetFillColor(kFillColor2); g->Draw("f");}
539  }
540  for (auto & g:g1) { g->SetFillColor(kWhite); g->DrawClone("f");}
541  for (auto & g:g1) { g->SetFillColor(kFillColor1); g->Draw("f");}
542  }
543  else{
544  //for (auto &g:g3) {
545  // g->SetFillColor(kFillColor3); g->Draw("f");}
546  //JoinGraphs(g3[0], g3[1], kWhite)->Draw("f");
547  JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f c");
548  JoinGraphs(g2[0], g2[1], kWhite)->Draw("f c");
549  JoinGraphs(g2[0], g2[1], kFillColor2)->Draw("f c");
550  /*for (auto &g:g2) {
551  g->SetFillColor(kWhite); g->DrawClone("f");
552  g->SetFillColor(kFillColor2); g->Draw("f");
553  }*/
554  JoinGraphs(g1[0], g1[1], kWhite)->Draw("f c");
555  JoinGraphs(g1[0], g1[1], kFillColor1)->Draw("f c");
556  }
557 
558  for (auto & gs:{g3,g2,g1}){
559  for (auto & g:gs) {
560  g->SetLineColor(kDarkColor);
561  g->SetLineWidth(3);
562  g->Draw("l");
563  }
564  }
565  }
566 
567 }
568 
569 void DrawReactorConstraint(double mean, double err)
570 {
571  TBox* box = new TBox(mean-err, 0, mean+err, 2);
572  box->SetFillColorAlpha(kGray,0.7);
573 
574  TLine* linLo = new TLine(mean-err, 0, mean-err, 2);
575  linLo->SetLineColor(kGray+2);
576  linLo->SetLineWidth(2);
577 
578  TLine* linHi = new TLine(mean+err, 0, mean+err, 2);
579  linHi->SetLineColor(kGray+2);
580  linHi->SetLineWidth(2);
581 
582  box->Draw();
583  linLo->Draw();
584  linHi->Draw();
585  gPad->RedrawAxis();
586 }
587 
588 // ssth23dmsq32 NH
589 // cafe joint_fit_2019_contours.C false true false joint_realData_both_onlyNH true false true true false true false
590 // ssth23dmsq32 IH
591 // cafe joint_fit_2019_contours.C false true false joint_realData_both_onlyIH true false true true false true false
592 
593 // deltassth23 NH
594 // cafe joint_fit_2019_contours.C false true false joint_realData_both_onlyNH false false true true false true false
595 // deltassth23 IH
596 // cafe joint_fit_2019_contours.C false true false joint_realData_both_onlyIH false false true true false true false
T max(const caf::Proxy< T > &a, T b)
TFile * MakeDataReleaseFileContour(TString preffix, bool dmsqSurf)
double minY
Definition: plot_hist.C:9
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
double GetValue(const osc::IOscCalcAdjustable *osc) const override
double th23
Definition: runWimpSim.h:98
double delta
Definition: runWimpSim.h:98
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Forward to wrapped Var&#39;s GetValue()
TH2 * Gaussian68Percent2D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 2D in gaussian approximation.
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 DrawHieTag(int hie, bool th13=true, bool vert=false)
void CylindricalSmoothing(TH2 *h)
TString hists[nhists]
Definition: bdt_com.C:3
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Definition: FitVars.cxx:42
osc::OscCalcDumb calc
void SetPaletteDeltaCyclicNew()
tuple white
Definition: rootlogon.py:71
string outfilename
knobs that need extra care
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)
string infile
virtual IOscCalc * Copy() const
Definition: OscCalcDumb.h:19
TH2 * Gaussian3Sigma2D(const FrequentistSurface &s)
Up-value surface for 3 sigma confidence in 2D in gaussian approximation.
void plot_joint_fit_2020_contours(bool corrSysts=false, bool runOnGrid=false, TString options="joint_fake2019_both_onlyIH", bool dmsqSurf=false, bool th13Surf=false, bool fccorr=false, bool zoomIn=true, bool fillContour=false, bool smooth=true, bool fc_overlay=false)
void GraphsToDataRelease(TFile *file, std::vector< TGraph * > graphs, TString label, TString hiestr)
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
void Add(const FCPoint &pt, std::string plot)
Definition: FCSurface.cxx:39
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
void BestFitToDataRelease(TFile *file, TMarker &bf, TString label, TString hiestr)
TLegend * ContourLegend(int hie, bool fillContour, bool fccorr, Int_t kFillColor1, Int_t kFillColor2, Int_t kFillColor3, Int_t kDarkColor, bool bestFit)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const Style_t k2SigmaConfidenceStyle
Definition: Style.h:71
double dmsq32
void DrawReactorConstraint(double mean, double err)
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
assert(nhit_max >=nhit_nbins)
const Float_t kBlessedLabelFontSize
Definition: Style.h:90
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
TH2 * SurfaceForSignificance(double sig)
Definition: FCSurface.cxx:103
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
static std::unique_ptr< FrequentistSurface > LoadFrom(TDirectory *dir, const std::string &name)
double th13
Definition: runWimpSim.h:98
static std::unique_ptr< FCCollection > LoadFromFile(const std::string &wildcard)
const std::string outdir
Pseudo-experiments, binned to match a log-likelihood surface.
Definition: FCSurface.h:14
surf
Definition: demo5.py:35
unsigned int uint