plot_joint_fit_2020_slices.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 "OscLib/IOscCalc.h"
18 #include "Utilities/rootlogon.C"
19 
20 #include "../joint_fit_2020_loader_tools.h"
21 #include "../joint_fit_2020_datarelease_tools.h"
22 
23 #include "TCanvas.h"
24 #include "TBox.h"
25 #include "TColor.h"
26 #include "TGraph.h"
27 #include "TVectorD.h"
28 #include "TF1.h"
29 #include "TLegend.h"
30 #include "TText.h"
31 #include "TLatex.h"
32 #include "TPad.h"
33 #include "TLine.h"
34 #include "TMarker.h"
35 #include "TStyle.h"
36 #include "TSystem.h"
37 #include "TGaxis.h"
38 #include "TGraphSmooth.h"
39 
40 #include <algorithm>
41 #include <vector>
42 #include <string>
43 #include <fstream>
44 #include <sstream>
45 
46 using namespace ana;
47 
48 
49 TGraph * signFuncGr;
50 Double_t signFunc(Double_t *x, Double_t *t);
51 void DrawSignMarker(double x, double y, int color);
52 void FindSignPoint(TGraph* gr, double sigma, double lowx, double highx, int color);
53 void HighlightSignPoints(TGraph* gr, TString grname);
54 void POTtag(bool both, bool FHCOnly, bool RHCOnly);
55 void AverageCyclicPoints(TGraph* g);
56 void AddCyclicPoints(TGraph* g);
57 TGraph* HistoToCyclicGraph(TGraph* g, bool useRoot);
58 void SmoothWithFourierFit(TH1* hist,
59  int nbins, double minX, double maxX,
60  int fourierFitN);
61 //why is this function so ugly
62 TGraph * SmoothWithFourierFit(TGraph* gr0, bool isDelta,
63  bool useRoot, int fourierFitN);
64 TGraph* FCCorrectSlice(TGraph* sqrtslice, TString fcFileName,
65  bool fourierFit, std::string plot_name,
66  bool useRoot, bool fccol);
67 
68 bool is_file_exist(const char *fileName)
69 {
70  std::ifstream infile(fileName);
71  return infile.good();
72 }
73 
74 
75 void plot_joint_fit_2020_slices(bool corrSysts = false,
76  bool runOnGrid = false,
77  TString options="joint_realData_both",
78  bool th23slice = false,
79  bool octantSlice = true,
80  bool dmsqSlice = false,
81  bool th13Slice = false,
82  bool fccorr = false,
83  bool fourierFit = true)
84 
85 {
86  bool nueOnly = options.Contains("nueOnly");
87  bool numuOnly = options.Contains("numuOnly");
88  bool joint = options.Contains("joint");
89  assert (nueOnly || numuOnly || joint);
90 
91  bool FHCOnly = options.Contains("FHCOnly");
92  bool RHCOnly = options.Contains("RHCOnly");
93  bool both = options.Contains("both");
94  assert (FHCOnly || RHCOnly || both);
95 
96  bool fake2019 = options.Contains("fake2019");
97  bool realData = options.Contains("realData");
98 
99  auto suffix = options;
100  if(corrSysts) suffix+="_systs";
101  else suffix+="_stats";
102 
103  if(th23slice) suffix+="_th23";
104  if(dmsqSlice) suffix+="_dmsq";
105  if(th13Slice) suffix+="_th13";
106  if(!octantSlice)suffix+="_noOct";
107 
108  TString outdir = "/nova/ana/3flavor/Ana2020/Fits/RealData/WithSeeds/";
109  if (runOnGrid) outdir = "./";
110 
111  TString outFCfilename (outdir + "slices_FCInput_2020_" + suffix);
112  if(joint && corrSysts && !th23slice && ! dmsqSlice && !th13Slice) outFCfilename = outdir + "slices_FCInput_2020_" +suffix;
113 
114  TString outfilename (outdir + "hist_slices_2020_" + suffix);
115 
116  TString infilename = outdir;
117 
118  infilename += "hist_slices_2020_";
119  infilename += options;
120  infilename += corrSysts?"_systs":"_stats";
121  if(th23slice) infilename+="_th23";
122  if(dmsqSlice) infilename+="_dmsq";
123  if(th13Slice) infilename+="_th13";
124  if(!octantSlice)infilename+="_noOct";
125  infilename += ".root";
126 
127  std::cout << "Making plots from " << infilename << std::endl;
128  TFile * infile = new TFile (infilename,"read");
129  TFile * official =0;
130  TString officialpreffix = "./";
131  bool makedatarelease = true;
132  if(makedatarelease &&
133  !(corrSysts && options.Contains("joint") && fccorr )) {
134  std::cerr << "\nWARNING: "
135  << "makedatarelease is set to true"
136  << "for a weird combination of options\n";
137  }
138  if(makedatarelease) {
139  official = MakeDataReleaseFileSlice(officialpreffix,
140  th23slice, dmsqSlice);
141  }
142  auto mins =* (TVectorD*)infile->Get("overall_minchi");
143  double minchi23 = mins[0];
144 
145  std::vector <TString> sliceNames = {};
146  if (!th23slice && ! dmsqSlice && !th13Slice) sliceNames = {"slice_delta_"};
147  if (th23slice) sliceNames = {"slice_th23_"};
148  if (dmsqSlice) sliceNames = {"slice_dmsq_"};
149  if (th13Slice) sliceNames = {"slice_th13_"};
150 
151  // it shouldn't play any role, but copy old style seeding here for sake of naming, everywhere is just .lable
152 
153  struct th23helper{
154  std::vector<double> seeds;
155  const IFitVar * var;
156  TString label;
157  };
158  std::vector <th23helper> th23seeds;
159  if(octantSlice) {
160  th23seeds.push_back( { {0.55}, &kFitSinSqTheta23UpperOctant, "UO"});
161  th23seeds.push_back( { {0.45}, &kFitSinSqTheta23LowerOctant, "LO"});
162  }
163  else th23seeds.push_back({ {0.45, 0.5, 0.55}, &kFitSinSqTheta23, ""});
164 
165  //end of old-style-seeding
166 
167  for(TString sliceName:sliceNames){
168  std::vector < std::pair <TGraph*, TString> > graphs;
169 
170  if(sliceName.Contains("delta")) {
171  for (TString hie:{"IH", "NH"}){
172  for(auto const & thseed:th23seeds){
173  //auto h = (TH1D*) infile ->Get( sliceName + hie + thseed.label); // previously was TH1 saved into the file
174  auto gr = (TGraph*) infile ->Get( sliceName + hie + thseed.label);
175  if(!gr) {
176  std::cerr << "Problem " << sliceName + hie + thseed.label << "\n";
177  exit(1);
178  }
179  if(fccorr)
180  graphs.push_back({gr, hie + " " + thseed.label});
181  else
182  graphs.push_back({HistoToCyclicGraph(gr, true), hie + " " + thseed.label});
183  }
184  }
185  DrawSliceCanvas(sliceName + (fccorr?"fccorr":""), 5, 0, 2);
186  }
187  else {
188  for (TString hie:{"NH","IH"}) {
189  for(auto const & thseed:th23seeds){
190  auto gr = (TGraph*) infile ->Get( sliceName + hie + thseed.label);
191  if(!gr) {std::cerr << "Problem " << sliceName + hie + thseed.label << "\n"; }
192  graphs.push_back({gr, hie + " " + thseed.label});
193  }
194  }
195 
196  if(sliceName.Contains("th23")) {
197  DrawSliceCanvas(sliceName + (fccorr?"fccorr":""), 3,
198  (nueOnly ? 0 : 0.32), (nueOnly ? 1 : 0.72));
199  }
200 
201  if(sliceName.Contains("dmsq")) {
202  DrawSliceCanvas(sliceName + (fccorr?"fccorr":""), 3, 2.15, 2.75);
203  }
204  if(sliceName.Contains("th13")) {
205  DrawSliceCanvas(sliceName + (fccorr?"fccorr":""), 3, 0.05, 0.2);
206  }
207  }
208 
209  TString fcFolder="/nova/ana/nu_e_ana/Ana2020/FC/slices/";
210  std::string plot_name = "";
211  if(th23slice) {
212  plot_name = "ssth23";
213  fcFolder += "th23/";
214  }
215  else if(th13Slice) {
216  plot_name = "th13";
217  fcFolder += "th13/";
218  }
219  else if(dmsqSlice) {
220  plot_name = "dmsq";
221  fcFolder += "dmsq/";
222  }
223  else {
224  plot_name = "delta";
225  fcFolder += "delta/";
226  }
227 
228  if(fccorr) {
229  if(!joint || !corrSysts)
230  { std::cerr << "\nWe don't have FC\n\n"; exit(1);}
231 
232  // TString fcFolder="./FC/";
233  for (auto &gr:graphs) {
234  TString fcSliceN = plot_name;
235  fcSliceN += gr.second.Contains("NH")? "_nh" : "_ih";
236  if(!th23slice)
237  fcSliceN += gr.second.Contains("UO")? "uo" : "lo";
238  fcSliceN += ".root";
239 
240  // only use fourier fit for period boundary conditions
241  // Fourier fit doesn't like delta_UO_nh. Would need
242  // lots of harmonics to fit it correctly anyway
243  if (is_file_exist(fcFolder+fcSliceN)){
244  std::cout<<"Found FC file for this curve"<<std::endl;
245  bool useFourier = (fcSliceN.Contains("delta") &&
246  !fcSliceN.Contains("delta_nhuo"));
247 
248  gr.first = FCCorrectSlice(gr.first, fcFolder + fcSliceN,
249  fourierFit, plot_name, !useFourier, true);
250  }
251  }
252 
253  }//end fccorr
254 
255  for (auto &gr:graphs) {
256  //if(gr.second.Contains("IH")) gr.first->SetLineColor(kPrimColorIH);
257  //if(gr.second.Contains("NH")) gr.first->SetLineColor(k3SigmaConfidenceColorNH);
258  if(gr.second.Contains("IH")) gr.first->SetLineColor(cih_dark);
259  if(gr.second.Contains("NH")) gr.first->SetLineColor(cnh);
260  if(gr.second.Contains("LO")) gr.first->SetLineStyle(7);
261  gr.first->SetLineWidth(3);
262  if(dmsqSlice){
263  for (int i =0; i < gr.first->GetN(); ++i){
264  gr.first->SetPoint(i,1000*fabs(gr.first->GetX()[i]),gr.first->GetY()[i]);
265  }
266  }
267  TGraphSmooth *test = new TGraphSmooth("normal");
268  TGraph * test_sm = test->SmoothKern(gr.first,"normal", 0.01);
269  test_sm->SetLineColor(gr.first->GetLineColor());
270  test_sm->SetLineStyle(gr.first->GetLineStyle());
271  test_sm->SetLineWidth(3);
272  test_sm->Draw("l same");
273  //gr.first->Draw("l same");
274  }
275 
276  POTtag(both, FHCOnly, RHCOnly);
277  PreliminarySide();
278  gPad->RedrawAxis();
279 
280  auto leg = SliceLegend(graphs, !dmsqSlice && !th23slice && !th13Slice, 0.7, dmsqSlice);
281  leg->Draw();
282 
283  // outdir = "test_delta_slice/";
284  //gPad->SetGrid();
285  for(TString ext: {".pdf",".root"}) {
286  gPad->Print("./slice_" + suffix +
287  // gPad->Print(TString("./test/") + [>(corrSysts?"syst/":"stat/")+<] "slice_" + suffix +
288  (fccorr?"_fccorr":"") +
289  (fourierFit?"_smooth":"") +
290  ext);
291  }//end ext
292 
293  bool annotate = 1;
294  if(annotate) {
295 
296  for (auto &gr:graphs) {
297  HighlightSignPoints(gr.first, sliceName);
298  }
299  gPad->SetGrid();
300  gPad->Print("./points_slice_" + suffix +
301  // gPad->Print(TString("./test/") + [>(corrSysts?"syst/":"stat/")+<] "points_slice_" + suffix +
302  (fccorr?"_fccorr":"") +
303  (fourierFit?"_smooth":"") +
304  ".pdf");
305  }
306  if(official) GraphsToDataRelease(official, graphs);
307  }//end slicename
308 
309 
310  bool debug = false;
311 
312  if(debug){
313 
314  TFile * inFCfile = new TFile(outFCfilename+".root","read");
315  std::cout << "\n\nOpen profiled " << outFCfilename+".root\n" << std::endl;
316 
317  std::vector <TString> strprof = {"DmSq32","SinSq2Theta13","SinSqTheta23"};
318  if(th23slice) strprof = {"DmSq32","SinSq2Theta13","DeltaCPpi"};
319  if(th13Slice) strprof = {"DmSq32","DeltaCPpi", "SinSqTheta23"};
320  if(dmsqSlice) strprof = {"SinSqTheta23","SinSq2Theta13","DeltaCPpi"};
321 
322  auto systs = GetJointFitSystematicList(corrSysts, nueOnly, numuOnly, FHCOnly||both, RHCOnly||both, true, true);;
323 
324  for (const ISyst* s : systs){
325  strprof.push_back(s->ShortName());
326  }
327 
328  for(auto const & str:strprof){
329  new TCanvas;
330  double miny=-1.5, maxy=1.5;
331  if(str=="DeltaCPpi") {miny=0; maxy=2;}
332  if(str=="DmSq32") {miny=2.2; maxy=2.9;}
333  if(str=="SinSq2Theta13") {miny=0.080; maxy=0.086;}
334  if(str=="SinSqTheta23") {miny=0.3; maxy=0.8;}
335 
336  auto axes = new TH2F("",";#delta_{CP}/#pi;" + str, 100, 0, 2, 100, miny, maxy);
337  if(dmsqSlice) axes = new TH2F("",";|#Deltam^{2}_{32}| (10^{-3} eV^{2});" + str, 100, 2.2, 2.8, 100, miny, maxy);
338  if(th23slice) axes = new TH2F("",";sin^{2}#theta_{23};" + str, 100, 0.3, 0.7, 100, miny, maxy);
339  if(th13Slice) axes = new TH2F("",";sin^{2}#theta_{13};" + str, 100, 0.01, 0.2, 100, miny, maxy);
340  axes->Draw("hist");
341 
342  for (TString hie:{"NH","IH"}){
343  for(auto const & thseed:th23seeds){
344  auto gr = (TGraph*) inFCfile->Get(hie + thseed.label + "_" + str);
345  if(!gr) {
346  std::cerr << "Problem " << hie + thseed.label + "_" + str << "\n";
347  exit(1);
348  }
349  if(str=="DmSq32"){
350  for (int i = 0; i < (int) gr->GetN(); ++i){
351  gr->SetPoint(i, gr->GetX()[i], 1000* fabs(gr->GetY()[i]));
352  }
353  }
354  if(dmsqSlice){
355  for (int i = 0; i < (int) gr->GetN(); ++i){
356  gr->SetPoint(i, 1000*fabs(gr->GetX()[i]), gr->GetY()[i]);
357  }
358  }
359  //gr->SetLineColor(hie=="NH"?k2SigmaConfidenceColorNH:kPrimColorIH);
360  gr->SetLineColor(hie=="NH"?cnh:cih_dark);
361  if(thseed.label.Contains("LO")) gr->SetLineStyle(7);
362  gr->SetLineWidth(3);
363  gr->SetMarkerColor(gr->GetLineColor());
364  gr->Draw("l same");
365  }//end th23
366  }//end hie
367  gPad->SetGrid();
368  gPad->Print("./debug_slice_prof_" + suffix + "_"+ str + ".pdf");
369  }//end profile vars
370  }//end debug
371 
372 }//end all
373 
374 
375 /// And supporting functions:
376 
377 Double_t signFunc(Double_t *x, Double_t *t){
378 
379  double y = signFuncGr->Eval(x[0]);
380  if(t[0] > 0 && abs(y - t[0]) > 0.05) return 999.;
381  return abs(y - t[0]);
382 }
383 void DrawSignMarker(double x, double y, int color)
384 {
385  TMarker * m = new TMarker (x, y,29);
386  m->SetMarkerColor(color);
387  m->SetMarkerSize(3);
388  TString pstr = TString::Format("#bf{(%.2f,%.2f)}",x,y);
389  TLatex * ltx = new TLatex(x,y,pstr.Data());
390  ltx->SetTextAngle(45);
391  ltx->SetTextColor(color);
392  ltx->Draw();
393  m->Draw();
394 }
395 void FindSignPoint(TGraph* gr, double sigma, double lowx=1, double highx=2, int color = kMagenta)
396 {
397  double minX = gr->GetX()[0];
398  double maxX = gr->GetX()[gr->GetN()];
399  signFuncGr = gr;
400  TF1 * f1 = new TF1("sign",signFunc,minX,maxX,1);
401  f1->SetParameter(0,sigma);
402  double xmin = f1->GetMinimumX(lowx,highx);
403 
404  if(f1->Eval(xmin) < 5) DrawSignMarker(xmin, gr->Eval(xmin), color);
405 }
406 
407 
408 
409 void HighlightSignPoints(TGraph* gr, TString grname)
410 {
411  struct SPoint{
412  double sigma;
413  double xlow;
414  double xhigh;
415  };
416  std::vector <SPoint > sig_range;
417  //fish for local minima points / 123 sgma points
418  if(grname.Contains("dmsq")){
419  sig_range = {{0.,2.3,2.6},
420  {1,2.3,2.4},{1,2.4,2.6},
421  {2,2.2,2.5},{2,2.5,2.7}};
422  }
423  if(grname.Contains("th23")){
424  sig_range = {{0.,0,0.5},{0,0.5,0.7},{0,0.5,0.5001},
425  {1,0.3,0.45},{1,0.45,0.5},{1,0.5,0.515},{1,0.515,0.55},{1,0.54,0.7},
426  {2,0.3,0.5},{2,0.5,0.7}};
427  }
428  if(grname.Contains("th13")){
429  sig_range = {{0.,0,0.1},{0,0.1,0.15},{0,0.15,0.2},
430  {1,0.0,0.1},{1,0.1,0.15},{1,0.15,0.2},
431  {2,0.0,0.1},{2,0.1,0.2}};
432  }
433  if(grname.Contains("delta")){
434  sig_range = {{0.,1., 2.}, {0.,0., 1.},
435  {1., 0, 0.5},{1.,0.5, 1.},{1.,1, 1.7},{1.,1.7, 2.},
436  {2., 0, 0.5},{2.,0.5, 1.},{2.,1, 1.3},{2.,1.5, 2.},
437  {3., 0, 0.5}, {3,0.5,1.5}, {3,1.7,2.}};
438 
439  }
440  int color = kMagenta;
441  for (auto &sr:sig_range){
442  if(sr.sigma == 1 ) color = kGreen + 1;
443  if(sr.sigma == 2 ) color = kAzure ;
444  if(sr.sigma == 3 ) color = kOrange ;
445  FindSignPoint(gr,sr.sigma,sr.xlow, sr.xhigh, color);
446  }
447 }
448 
449 void POTtag(bool both, bool FHCOnly, bool RHCOnly)
450 {
451  TLatex* ltx = new TLatex();
452  auto ltx1 = new TLatex(.15, 0.93, TString::Format("NOvA FD"));
453  if (FHCOnly){ ltx = new TLatex(.42, 0.93, TString::Format("%g#times10^{20} POT equiv. of #nu data only", kAna2020FHCPOT/1E20).Data());}
454  if (RHCOnly){ ltx = new TLatex(.52, 0.93, TString::Format("%2.4g#times10^{20} POT of #bar{#nu} data only", kAna2020RHCPOT/1E20).Data());}
455  if(both) {ltx = new TLatex(.32, 0.93, TString::Format("%g#times10^{20} POT equiv #nu + %2.4g#times10^{20} POT #bar{#nu}", kAna2020FHCFullDetEquivPOT/1E20, kAna2020RHCPOT/double(1E20)).Data());}
456  ltx1->SetNDC();
457  ltx1->SetTextAlign(11);
458  ltx1->SetTextSize(0.8*kBlessedLabelFontSize);
459  ltx1->Draw();
460  ltx->SetNDC();
461  ltx->SetTextAlign(11);
462  ltx->SetTextSize(0.8*kBlessedLabelFontSize);
463  ltx->Draw("same");
464 }
465 
466 void AverageCyclicPoints(TGraph* g)
467 {
468  double x0, y0;
469  double x2pi, y2pi;
470  g->GetPoint(g->GetN()-1, x2pi, y2pi);
471  g->GetPoint(1, x0, y0);
472  std::cout << "(x0, y0) = (" << x0 << ", " << y0 << ")\n";
473  std::cout << "(x2pi, y2pi) = (" << x2pi << ", " << y2pi << ")\n";
474 
475  double avgY = ( y0 + y2pi ) / 2;
476  std::cout << "Average = " << avgY << std::endl;
477  g->SetPoint(0, 0, avgY);
478  g->SetPoint(g->GetN(), 2, avgY);
479 
480  std::cout << "\nFound average cyclic point: (0, " << avgY << ")\n";
481  std::cout << "\nAdded average cyclic points " << 0 << ": " << g->GetX()[0] << " " << g->GetY()[0];
482  std::cout << ", " << g->GetN() << ": " << g->GetX()[g->GetN()-1] << " " << g->GetY()[g->GetN()-1] << "\n\n";
483 }
484 
485 
486 
487 void AddCyclicPoints(TGraph* g)
488 {
489  double x, y;
490  g->GetPoint(g->GetN()-1, x, y);
491  g->SetPoint(0, x-2, y);
492  g->GetPoint(1, x, y);
493  g->SetPoint(g->GetN(), x+2, y);
494 
495  std::cout << "\n Added cyclic points " << 0 << " " << g->GetX()[0]<< " " << g->GetY()[0]
496  << ", " << g->GetN()-1 << x+2 << y << "\n\n";
497 }
498 
499 
500 TGraph* HistoToCyclicGraph(TGraph* g, bool useRoot = false)
501 {
502  // TGraph* g = new TGraph;
503  // g->SetPoint(0, 0, 0); // temporary
504  // for(int i = 1; i <= h->GetNbinsX(); ++i)
505  // g->SetPoint(i, h->GetXaxis()->GetBinCenter(i), h->GetBinContent(i));
506 
507  //g->SetLineColor(h->GetLineColor());
508  g->SetLineWidth(2);
509 
510  if(useRoot)
512  else
513  AddCyclicPoints(g);
514 
515  return g;
516 }
517 
518 
520  int nbins = 80, double minX = 0, double maxX = 2,
521  int fourierFitN=3)
522 {
523  FitToFourier fitF (hist, minX, maxX, fourierFitN);
524  auto f1= fitF.Fit();
525  hist->Reset();
526  hist->Add(f1);
527 }
528 //why is this function so ugly
529 TGraph * SmoothWithFourierFit(TGraph* gr0,
530  bool isDelta = true,
531  bool useRoot = false,
532  int fourierFitN = 2)
533 {
534  int nbins = gr0->GetN();
535  double minX = gr0->GetX()[0];
536  double maxX = gr0->GetX()[nbins-1];
537  //TGraph to center bin
538  minX = minX - (gr0->GetX()[1] - gr0->GetX()[0])/2;
539  maxX = maxX + (gr0->GetX()[1] - gr0->GetX()[0])/2;
540 
541  //Fourier fit relies on periodicity, ignore extra points
542  if(isDelta && !useRoot){
543  nbins = 60;
544  minX = 0;
545  maxX = 2;
546  }
547  std::cout << "Smooth " << nbins << " nbins " << minX << ", " << maxX << "\n";
548 
549  TH1D * hist = new TH1D(UniqueName().c_str(),"",nbins,minX,maxX);
550 
551  for (int i = 0; i < nbins; ++i){
552  double val = gr0->Eval(hist->GetBinCenter(i+1));
553  hist->SetBinContent(i + 1, val * val); //use dchisq not chi2 to smooth
554  }
555  if(useRoot) hist->Smooth(2);
556  else SmoothWithFourierFit(hist, nbins, minX, maxX, fourierFitN);
557  //now return sqrt
558  for (int i = 1; i <= nbins; ++i){
559  hist->SetBinContent(i, sqrt(hist->GetBinContent(i)));
560  }
561  /*if(isDelta && !useRoot) // NB: commented, previously HistoToCyclicGraph used TH1D, now there is a graph
562  return HistoToCyclicGraph(hist, useRoot); // some corrections are needed
563  else */return new TGraph(hist);
564 }
565 
566 TGraph* FCCorrectSlice(TGraph* sqrtslice, TString fcFileName,
567  bool fourierFit, std::string plot_name,
568  bool useRoot = false, bool fccol=false)
569 {
570  bool th23plot = fcFileName.Contains("th23");
571  bool dmsqplot = fcFileName.Contains("dmsq");
572  bool dmsqnhuo = fcFileName.Contains("dmsq_nhuo");
573  TString hiername = fcFileName.Contains("nh")? "NH":"IH";
574 
575  double minX = 0;
576  double maxX = 2;
577  if (th23plot) {minX = 0.3; maxX = 0.7;}
578  if (dmsqplot) {
579  minX = hiername=="NH" ? 2e-3 : -3e-3;
580  maxX = hiername=="NH" ? 3e-3 : -2e-3;
581  }
582 
583  //We want slice, not sqrt
584  auto n= sqrtslice->GetN();
585 
586  auto sliceXH = new TGraph(n);
587  for (int i = 0; i < n;i++ ){
588  sliceXH->SetPoint(i,sqrtslice->GetX()[i],util::sqr(sqrtslice->GetY()[i]));
589  }
590  FCSurface *fcXH = 0;
591  if(!fccol) fcXH = FCSurface::LoadFromFile(fcFileName.Data()).release();
592  else {
593  auto fccol = FCCollection::LoadFromFile(fcFileName.Data()).release();
594  //Only one Y bin required for slices.
595  // All slices this year had 60 bins. Need to find a better way
596  // next year
597  fcXH = new FCSurface(60,minX,maxX,1,0,1);
598  fcXH->Add(*fccol, plot_name);
599  }
600 
601 
602  TGraph* sigXH = new TGraph;
603  TGraph* pcXH = new TGraph;
604 
605  TFile* bestfitfile = new TFile("/nova/ana/nu_e_ana/Ana2019/Results/BestFits/bestfits_joint_realData_both_systs.root", "read");
606  auto calcbf = (osc::IOscCalcAdjustable*)(LoadFrom<osc::IOscCalc>(bestfitfile, "NH_UO_calc")).release();
607  double dmsq32bf = kFitDmSq32.GetValue(calcbf);
608  double th23bf = kFitSinSqTheta23.GetValue(calcbf);
609  double deltabf = kFitDeltaInPiUnits.GetValue(calcbf);
610  for(int i = 0; i < n; ++i){
611  const double delta = sliceXH->GetX()[i];
612  const double upXH = sliceXH->GetY()[i];
613 
614  // Fall back to gaussian up if dchisq is too large or FCBin is empty (mostly for deltaCP IH slices)
615  if(upXH > 13 || fcXH->GetBin(i+1,1)->Empty()){
616  std::cout << " Falling back on gaussian up value: " << upXH
617  << " for bin " << i << " , delta = " << delta << std::endl;
618  sigXH->SetPoint(sigXH->GetN(), delta, sqrt(upXH));
619  }
620  else {
621  double pXH;
622  if(upXH > 0) pXH = fcXH->GetBin(i+1, 1)->SignificanceForUpValue(upXH);
623  else pXH = 0; //HACK HACK
624  sigXH->SetPoint(sigXH->GetN(), delta, PValueToSigma(1-pXH));
625  }
626  if(dmsqnhuo && (i+1 < n) && (delta < dmsq32bf) && (dmsq32bf < sliceXH->GetX()[i+1])){
627  sigXH->SetPoint(sigXH->GetN(), dmsq32bf, 0.);
628  }
629  if ( special_smoothing && (i+1 < n) && (delta < deltabf) && (deltabf < sliceXH->GetX()[i+1])){
630  std::cout << "Setting BFP at deltaCP = " << deltabf << std::endl;
631  sigXH->SetPoint(sigXH->GetN(), deltabf, 0.);
632  }
633  if ( fcFileName.Contains("th23_nh") && (i+1 < n) && (delta < th23bf) && (th23bf < sliceXH->GetX()[i+1])){
634  std::cout << "Setting BFP at th23 = " << th23bf << std::endl;
635  sigXH->SetPoint(sigXH->GetN(), th23bf, 0.);
636  }
637  }
638  if(fourierFit && !th23plot) {
639  sigXH = SmoothWithFourierFit(sigXH, fcFileName.Contains("delta"), useRoot);
640  }
641 
642  return sigXH;
643 }
644 
645 // th23
646 // cafe joint_fit_2019_slices.C false true false joint_realData_both true false false false true true
647 
648 // dmsq
649 // cafe joint_fit_2019_slices.C false true false joint_realData_both false true true false true true
650 
651 // delta
652 // cafe joint_fit_2019_slices.C false true false joint_realData_both false true false false true true
void DrawSliceCanvas(TString slicename, const double ymax, const double xmin=0, const double xmax=2.)
enum BeamMode kOrange
fileName
Definition: plotROC.py:78
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void HighlightSignPoints(TGraph *gr, TString grname)
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
double GetValue(const osc::IOscCalcAdjustable *osc) const override
double delta
Definition: runWimpSim.h:98
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Forward to wrapped Var&#39;s GetValue()
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:18
double maxy
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Definition: FitVars.cxx:171
T sqrt(T number)
Definition: d0nt_math.hpp:156
const FitSinSqTheta23UpperOctant kFitSinSqTheta23UpperOctant
Definition: FitVars.cxx:16
OStream cerr
Definition: OStream.cxx:7
void abs(TH1 *hist)
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
const FCBin * GetBin(int x, int y) const
Definition: FCSurface.cxx:152
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
const double kAna2020FHCFullDetEquivPOT
Definition: Exposures.h:234
double SignificanceForUpValue(double up) const
Definition: FCBin.cxx:59
const char * label
static std::unique_ptr< FCSurface > LoadFromFile(const std::string &fname)
Definition: FCSurface.cxx:170
double PValueToSigma(double p)
Compute the equivalent number of gaussian sigma for a p-value.
const XML_Char * s
Definition: expat.h:262
string outfilename
knobs that need extra care
const int nbins
Definition: cellShifts.C:15
string infile
caf::StandardRecord * sr
Float_t f1
bool is_file_exist(const char *fileName)
void GraphsToDataRelease(TFile *file, std::vector< TGraph * > graphs, TString label, TString hiestr)
TLegend * SliceLegend(std::vector< std::pair< TGraph *, TString > > &graphs, bool isDelta)
void Add(const FCPoint &pt, std::string plot)
Definition: FCSurface.cxx:39
const double kAna2020FHCPOT
Definition: Exposures.h:233
void AverageCyclicPoints(TGraph *g)
void DrawSignMarker(double x, double y, int color)
const double kAna2020RHCPOT
Definition: Exposures.h:235
double sigma(TH1F *hist, double percentile)
TF1 * Fit() const
Definition: Utilities.cxx:427
Double_t signFunc(Double_t *x, Double_t *t)
And supporting functions:
void AddCyclicPoints(TGraph *g)
OStream cout
Definition: OStream.cxx:6
TFile * MakeDataReleaseFileSlice(TString preffix, bool th23Slice, bool dmsqSlice)
const FitSinSqTheta23LowerOctant kFitSinSqTheta23LowerOctant
Definition: FitVars.cxx:17
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
TGraph * FCCorrectSlice(TGraph *sqrtslice, TString fcFileName, bool fourierFit, std::string plot_name, bool useRoot, bool fccol)
void SmoothWithFourierFit(TH1 *hist, int nbins, double minX, double maxX, int fourierFitN)
exit(0)
assert(nhit_max >=nhit_nbins)
void FindSignPoint(TGraph *gr, double sigma, double lowx, double highx, int color)
Interface definition for fittable variables.
Definition: IFitVar.h:16
void plot_joint_fit_2020_slices(bool corrSysts=false, bool runOnGrid=false, TString options="joint_realData_both", bool th23slice=false, bool octantSlice=true, bool dmsqSlice=false, bool th13Slice=false, bool fccorr=false, bool fourierFit=true)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
const Float_t kBlessedLabelFontSize
Definition: Style.h:90
void PreliminarySide()
Definition: rootlogon.C:114
Float_t e
Definition: plot.C:35
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
enum BeamMode kGreen
bool Empty() const
Definition: FCBin.h:23
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)
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:29
void POTtag(bool both, bool FHCOnly, bool RHCOnly)
const std::string outdir
TGraph * HistoToCyclicGraph(TGraph *g, bool useRoot)
double GetX(int ndb=14, int db=1, int feb=0, int pix=0)
Definition: PlotOnMon.C:111
Pseudo-experiments, binned to match a log-likelihood surface.
Definition: FCSurface.h:14
TGraph * signFuncGr
enum BeamMode string