joint_fit_2017_slices.C
Go to the documentation of this file.
3 #include "CAFAna/Fit/Fit.h"
12 #include "CAFAna/FC/FCSurface.h"
16 #include "CAFAna/Vars/FitVars.h"
17 #include "OscLib/IOscCalc.h"
19 
20 #include "TCanvas.h"
21 #include "TBox.h"
22 #include "TColor.h"
23 #include "TGraph.h"
24 #include "TVectorD.h"
25 #include "TF1.h"
26 #include "TLegend.h"
27 #include "TLatex.h"
28 #include "TLine.h"
29 #include "TMarker.h"
30 #include "TStyle.h"
31 #include "TSystem.h"
32 #include "TGaxis.h"
33 
34 #include <algorithm>
35 #include <vector>
36 #include <string>
37 
38 using namespace ana;
39 
40 
41 TGraph * signFuncGr;
42 Double_t signFunc(Double_t *x, Double_t *t){
43 
44  double y = signFuncGr->Eval(x[0]);
45  if(t[0] > 0 && abs(y - t[0]) > 0.05) return 999.;
46  return abs(y - t[0]);
47 }
48 void DrawSignMarker(double x, double y, int color)
49 {
50  TMarker * m = new TMarker (x, y,29);
51  m->SetMarkerColor(color);
52  m->SetMarkerSize(3);
53  TString pstr = TString::Format("#bf{(%.2f,%.2f)}",x,y);
54  TLatex * ltx = new TLatex(x,y,pstr.Data());
55  ltx->SetTextAngle(45);
56  ltx->SetTextColor(color);
57  ltx->Draw();
58  m->Draw();
59 }
60 void FindSignPoint(TGraph* gr, double sigma, double lowx=1, double highx=2, int color = kMagenta)
61 {
62  double minX = gr->GetX()[0];
63  double maxX = gr->GetX()[gr->GetN()];
64  signFuncGr = gr;
65  TF1 * f1 = new TF1("sign",signFunc,minX,maxX,1);
66  f1->SetParameter(0,sigma);
67  double xmin = f1->GetMinimumX(lowx,highx);
68 
69  if(f1->Eval(xmin) < 5) DrawSignMarker(xmin, gr->Eval(xmin), color);
70 }
71 
72 void HighlightSignPoints(TGraph* gr, TString grname)
73 {
74  struct SPoint{
75  double sigma;
76  double xlow;
77  double xhigh;
78  };
79  std::vector <SPoint > sig_range;
80  //fish for local minima points / 123 sgma points
81  if(grname.Contains("dmsq")){
82  sig_range = {{0.,2.3,2.6},
83  {1,2.3,2.4},{1,2.4,2.6},
84  {2,2.2,2.5},{2,2.5,2.7}};
85  }
86  if(grname.Contains("th23")){
87  sig_range = {{0.,0,0.5},{0,0.5,0.7},{0,0.5,0.5001},
88  {1,0.3,0.45},{1,0.45,0.5},{1,0.5,0.515},{1,0.515,0.55},{1,0.55,0.7},
89  {2,0.3,0.5},{2,0.5,0.7}};
90  }
91  if(grname.Contains("delta")){
92  sig_range = {{0.,1., 2.},
93  {1., 0, 0.5},{1.,0.5, 1.},{1.,1, 1.7},{1.,1.7, 2.},
94  {2., 0, 0.5},{2.,0.5, 1.},{2.,1, 1.3},{2.,1.5, 2.},
95  {3., 0, 0.5}, {3,0.5,1.5}, {3,1.7,2.}};
96 
97  }
98  int color = kMagenta;
99  for (auto &sr:sig_range){
100  if(sr.sigma == 1 ) color = kGreen + 1;
101  if(sr.sigma == 2 ) color = kAzure ;
102  if(sr.sigma == 3 ) color = kOrange ;
103  FindSignPoint(gr,sr.sigma,sr.xlow, sr.xhigh, color);
104  }
105 }
106 
107 
108 void MakeMaps ( std::vector <const IFitVar*> profvars,
109  std::map<const IFitVar*, TGraph*> &profMap,
110  std::vector <const ISyst* > systs,
111  std::map<const ISyst*, TGraph*> &systMap)
112 {
113 
114  for (const IFitVar * var:profvars)
115  profMap.insert(std::pair<const IFitVar*, TGraph*> (var, new TGraph()));
116  for (const ISyst* syst : systs)
117  systMap.insert(std::pair<const ISyst*, TGraph*> (syst, new TGraph()));
118 }
119 
120 void SaveMaps ( TDirectory * dir,
121  TString prefix,
122  std::vector <const IFitVar*> profvars,
123  std::vector <TString > profnames,
124  std::map<const IFitVar*, TGraph*> &profMap,
125  std::vector <const ISyst* > systs,
126  std::map<const ISyst*, TGraph*> &systMap)
127 {
128  TDirectory *tmp = gDirectory;
129  dir->cd();
130  for (int i = 0; i < (int) profvars.size(); ++i){
131  profMap[profvars[i]]->Write((prefix + "_" + profnames[i]));
132  }
133  for (const ISyst* s : systs)
134  systMap[s]->Write((prefix + "_" + s->ShortName()));
135  tmp->cd();
136 }
137 
138 
139 void DrawSliceCanvas(TString slicename, const double ymax,
140  const double xmin = 0, const double xmax = 2.){
141  auto c1 = new TCanvas(ana::UniqueName().c_str());
142  c1->SetFillStyle(0);
143  c1->SetLeftMargin(0.14);
144  c1->SetBottomMargin(0.15);
145  TString title;
146  if(slicename.Contains("delta")) {
147  title = ";delta;Significance (#sqrt{#Delta#chi^{2}})";
148  }
149  if(slicename.Contains("th23")) {
150  title = ";sin^{2}#theta_{23};Significance (#sqrt{#Delta#chi^{2}})";
151  }
152  if(slicename.Contains("dmsq")) {
153  title = ";|#Deltam^{2}_{32}| (10^{-3} eV^{2});Significance (#sqrt{#Delta#chi^{2}})";
154  }
155  auto axes = new TH2F("",title, 100, xmin, xmax, 100, 0, ymax);
156  if(slicename.Contains("fccorr")) axes->GetYaxis()->SetTitle("Significance (#sigma)");
158  axes->Draw();
159  axes->GetXaxis()->SetTitleSize(kBlessedTitleFontSize);
160  axes->GetYaxis()->SetTitleSize(kBlessedTitleFontSize);
161  axes->GetXaxis()->SetLabelSize(kBlessedLabelFontSize);
162  axes->GetYaxis()->SetLabelSize(kBlessedLabelFontSize);
163  axes->GetXaxis()->SetTitleOffset(0.8);
164  axes->GetYaxis()->SetTitleOffset(0.75);
165  if(slicename.Contains("delta")) XAxisDeltaCPLabels(axes);
166  c1->RedrawAxis();
167 }
168 
169 TLegend * SliceLegend(std::vector <std::pair <TGraph*, TString > > & graphs, bool isDelta)
170 {
171  TLegend * leg;
172  if(isDelta) leg = new TLegend(0.7, 0.6, 0.89,0.89);
173  else leg = new TLegend(0.7, 0.2, 0.89, 0.49);
174  leg->SetTextSize(kBlessedLabelFontSize);
175  leg->SetFillStyle(0);
176  for (auto &gr:graphs){
177  leg->AddEntry(gr.first,gr.second, "l");
178  }
179  return leg;
180 }
181 
182 
183 void POTtag()
184 {
185  TLatex* ltx = new TLatex(.16, .81, TString::Format("#splitline{NOvA FD}{%g#times10^{20} POT equiv.}", kAna2017FullDetEquivPOT/1E20).Data());
186  ltx->SetNDC();
187  ltx->SetTextAlign(11);
188  ltx->SetTextSize(kBlessedLabelFontSize);
189  ltx->Draw();
190 }
191 
192 void AddCyclicPoints(TGraph* g)
193 {
194  double x, y;
195  g->GetPoint(g->GetN()-1, x, y);
196  g->SetPoint(0, x-2, y);
197  g->GetPoint(1, x, y);
198  g->SetPoint(g->GetN(), x+2, y);
199 
200  std::cout << "\n Added cyclic points " << 0 << " " << g->GetX()[0]<< " " << g->GetY()[0]
201  << ", " << g->GetN() << x+2 << y << "\n\n";
202 }
203 
204 
205 TGraph* HistoToCyclicGraph(TH1* h)
206 {
207  TGraph* g = new TGraph;
208  g->SetPoint(0, 0, 0); // temporary
209  for(int i = 1; i <= h->GetNbinsX(); ++i)
210  g->SetPoint(i, h->GetXaxis()->GetBinCenter(i), h->GetBinContent(i));
211 
212  g->SetLineColor(h->GetLineColor());
213  g->SetLineWidth(2);
214 
215  AddCyclicPoints(g);
216 
217  return g;
218 }
220  int nbins = 80, double minX = 0, double maxX = 2,
221  int fourierFitN=3)
222 {
223  FitToFourier fitF (hist, minX, maxX, fourierFitN);
224  auto f1= fitF.Fit();
225  hist->Reset();
226  hist->Add(f1);
227 }
228 //why is this function so ugly
229 TGraph * SmoothWithFourierFit(TGraph* gr0,
230  bool isDelta = true,
231  bool useRoot = false,
232  int fourierFitN = 2)
233 {
234  int nbins = gr0->GetN();
235  double minX = gr0->GetX()[0];
236  double maxX = gr0->GetX()[nbins-1];
237  //TGraph to center bin
238  minX = minX - (gr0->GetX()[1] - gr0->GetX()[0])/2;
239  maxX = maxX + (gr0->GetX()[1] - gr0->GetX()[0])/2;
240 
241  //Fourier fit relies on periodicity, ignore extra points
242  if(isDelta && !useRoot){
243  nbins = 80;
244  minX = 0;
245  maxX = 2;
246  }
247  std::cout << "Smooth " << nbins << " nbins " << minX << ", " << maxX << "\n";
248 
249  TH1D * hist = new TH1D(UniqueName().c_str(),"",nbins,minX,maxX);
250 
251  for (int i = 0; i < nbins; ++i){
252  double val = gr0->Eval(hist->GetBinCenter(i+1));
253  hist->SetBinContent(i + 1, val * val); //use dchisq not chi2 to smooth
254  }
255  if(useRoot) hist->Smooth(2);
256  else SmoothWithFourierFit(hist, nbins, minX, maxX, fourierFitN);
257  // hist->Smooth(10);
258 
259  //now return sqrt
260  for (int i = 1; i <= nbins; ++i){
261  hist->SetBinContent(i, sqrt(hist->GetBinContent(i)));
262  }
263  if(isDelta && !useRoot) return HistoToCyclicGraph(hist);
264  else return new TGraph(hist);
265 }
266 
267 TGraph* FCCorrectSlice(TGraph* sqrtslice, TString fcFileName,
268  bool fourierFit, bool useRoot = false, bool fccol=false)
269 {
270  bool th23plot = fcFileName.Contains("th23");
271  bool dmsqplot = fcFileName.Contains("dmsq");
272  TString hiername = fcFileName.Contains("NH")? "NH":"IH";
273 
274  double minX = 0;
275  double maxX = 2;
276  if (th23plot) {minX = 0.3; maxX = 0.7;}
277  if (dmsqplot) {
278  minX = hiername=="NH" ? 2e-3 : -3e-3;
279  maxX = hiername=="NH" ? 3e-3 : -2e-3;
280  }
281 
282  //We want slice, not sqrt
283  auto n= sqrtslice->GetN();
284  auto sliceXH = new TGraph(n);
285  for (int i = 0; i < n;i++ ){
286  sliceXH->SetPoint(i,sqrtslice->GetX()[i],util::sqr(sqrtslice->GetY()[i]));
287  }
288  FCSurface *fcXH = 0;
289  if(!fccol) fcXH = FCSurface::LoadFromFile(fcFileName.Data()).release();
290  else {
291  auto fccol = FCCollection::LoadFromFile(fcFileName.Data()).release();
292  fcXH = new FCSurface(n,minX,maxX,3,0,1);
293  fcXH->Add(*fccol);
294  }
295  TGraph* sigXH = new TGraph;
296  TGraph* pcXH = new TGraph;
297 
298  for(int i = 0; i < n; ++i){
299  const double delta = sliceXH->GetX()[i];
300  const double upXH = sliceXH->GetY()[i];
301 
302  //HACK HACK
303  if((!fcFileName.Contains("delta") && upXH > 11)){
304  sigXH->SetPoint(sigXH->GetN(), delta, sqrt(upXH));
305  }
306  else {
307  double pXH;
308  if(upXH > 0) pXH = fcXH->GetBin(i+1, 2)->SignificanceForUpValue(upXH);
309  else pXH = 0; //HACK HACK
310  sigXH->SetPoint(sigXH->GetN(), delta, PValueToSigma(1-pXH));
311  }
312  }
313  if(fourierFit) {
314  sigXH = SmoothWithFourierFit(sigXH, fcFileName.Contains("delta"), useRoot);
315  }
316 
317  return sigXH;
318 }
319 
321  bool createFile=false,
322  bool corrSysts = false,
323  TString options="joint",
324  bool th23slice = false,
325  bool octantSlice = true,
326  bool dmsqSlice = false,
327  bool fccorr=false,
328  bool fourierFit = false)
329 {
330  bool nueOnly = options.Contains("nueOnly");
331  bool numuOnly = options.Contains("numuOnly");
332  bool joint = options.Contains("joint");
333  assert (nueOnly || numuOnly || joint);
334 
335  bool fakeNHLO = options.Contains("fakeNHLO");
336  bool fakeNHUO = options.Contains("fakeNHUO");
337  bool fake2017 = options.Contains("fake2017");
338  bool realData = options.Contains("realData");
339 
340  auto suffix = options;
341  suffix += decomp;
342  if(corrSysts) suffix+="_systs";
343  else suffix+="_stats";
344 
345  if(th23slice) suffix+="_th23";
346  if(dmsqSlice) suffix+="_dmsq";
347  if(!octantSlice)suffix+="_noOct";
348 
349  TString outdir = "/nova/ana/nu_e_ana/Ana2017/Results/slices/";
350  // TString outdir = "slices/";
351  TString outfilename (outdir + "hist_slices_2017_" + suffix);
352  TString outFCfilename (outdir + "slices_FCInput_2017_" + suffix);
353  if(joint && corrSysts && !th23slice && ! dmsqSlice)
354  outFCfilename = outdir + "slices_FCInput_2017_" + decomp;
355 
356  struct th23helper{
357  std::vector<double> seeds;
358  const IFitVar * var;
359  TString label;
360  };
361  std::vector <th23helper> th23seeds;
362  th23helper temp;
363  temp.seeds = {0.3};
364  temp.var = kFitSinSqTheta23LowerOctant;
365  temp.label = "LO";
366  if(octantSlice) {
367  th23seeds.push_back( { {0.45}, kFitSinSqTheta23LowerOctant, "LO"});
368  th23seeds.push_back( { {0.55}, kFitSinSqTheta23UpperOctant, "UO"});
369  }
370  else th23seeds.push_back({ {0.45,0.55}, &kFitSinSqTheta23, ""});
371 
372  if(createFile){
373 
374  //////////////////////////////////////////////////
375  // Load Nue and Numu experiments
376  //////////////////////////////////////////////////
377 
378  std::vector <const IPrediction * > preds;
379  std::vector <std::pair <TH1D*, double > > cosmics;
380  std::vector <Spectrum * > data;
381  std::vector <const IExperiment * > expts;
382 
383  auto calc_fake = DefaultOscCalc();
384  if(fakeNHLO) SetFakeCalc(calc_fake, 0.404, 2.7e-3, 1.48);
385  else if(fakeNHUO) SetFakeCalc(calc_fake, 0.623, 2.7e-3, 0.74);
386  else if(fake2017) SetFakeCalc(calc_fake, 0.56, 2.43e-3, 1.21);
387  else if(!realData) {std::cerr << "need setting for data\n"; exit(1);}
388 
389  if(!numuOnly) {
390  preds.push_back(GetNuePrediction2017(decomp, DefaultOscCalc(), corrSysts));
391  cosmics.push_back(GetNueCosmics2017());
392  }
393 
394  int nnumu = 4;
395 
396  if(!nueOnly) {
397  auto numu_preds = GetNumuPredictions2017(nnumu, corrSysts);
398  preds.insert(preds.end(),numu_preds.begin(), numu_preds.end());
399  auto numu_cosmics = GetNumuCosmics2017(nnumu);
400  cosmics.insert(cosmics.end(),numu_cosmics.begin(), numu_cosmics.end());
401  }
402 
403  if(realData){
404  if(!numuOnly){
405  data.push_back(GetNueData2017());
406  }
407  if(!nueOnly){
408  auto numu_data = GetNumuData2017(nnumu);
409  data.insert(data.end(),numu_data.begin(), numu_data.end());
410  }
411  }
412  for(int i = 0; i < (int) preds.size(); ++i){
413  auto thisdata = GetFakeData (preds[i],calc_fake, kAna2017POT, cosmics[i].first);
414  if(realData) thisdata = data[i];
415  expts.push_back(new SingleSampleExperiment(preds[i],*thisdata,
416  cosmics[i].first,cosmics[i].second));
417  }
418 
419  ////////////////////////////////////////////////////////////
420  // Add constraints, make experiments
421  ////////////////////////////////////////////////////////////
422  std::cout << "\nCreating multiexperiment\n" << std::endl;
423 
424  std::cout << "Adding WorldReactorConstraint2017\n";
425  expts.push_back(WorldReactorConstraint2017());
426  if(nueOnly) {
427  std::cout << "Adding Dmsq32ConstraintPDG2017\n";
428  expts.push_back(&kDmsq32ConstraintPDG2017);
429  }
430  std::cout << "Creating Multiexperiment with a total of "
431  << expts.size() << " experiments\n\n" << std::endl;
432  auto exptThis = new MultiExperiment(expts);
433 
434  ////////////////////////////////////////////////////////////
435  // Systematics
436  ////////////////////////////////////////////////////////////
437  std::cout << "Systematics for the fit:\n";
438  auto systs = GetJointFitSystematicList(corrSysts, !numuOnly, !nueOnly, true);
439 
440  std::cout << "\n\nSystematic correlations...\n";
441  if(!nueOnly && ! numuOnly && corrSysts){
442  exptThis->SetSystCorrelations(0, GetCorrelations(true));
443  auto notfornumu = GetCorrelations(false);
444  for(int i =0; i < nnumu; ++i) exptThis->SetSystCorrelations(i+1, notfornumu);
445  }
446 
447  ////////////////////////////////////////////////////////////
448  // Fit
449  ////////////////////////////////////////////////////////////
450  std::cout << "Starting the fit" << std::endl;
451 
453 
454  SystShifts auxShifts = SystShifts::Nominal();
455 
456  // In case some systs need different seeds
457  std::vector <SystShifts> seedShifts = {};
458  if(corrSysts){
459  for (double systshift:{+0.5,-0.5}){
460  SystShifts tempShifts (kMECq0ShapeSyst2017,systshift);
461  seedShifts.push_back(tempShifts);
462  }
463  }
464 
465  // Find the best fit points
466  double minchi23 = 1E20;
467 
468  for (auto & thseed: th23seeds){
469  for(int hie:{-1, 1}){
470 
471  std::cout << "\n\nFinding best fit " << (hie > 0 ? "NH " : "IH ")
472  << thseed.label << ", "
473  << "ssth23 seeds ";
474  for (auto const & s:thseed.seeds) std::cout << s << ", ";
475  std::cout << std::endl;
476 
477  std::vector <const IFitVar*> fitvars = {&kFitDeltaInPiUnits,
478  thseed.var,
481 
482  MinuitFitter fit23(exptThis, fitvars, systs);
483 
484  auto thisminchi = fit23.Fit(calc, auxShifts,
485  {
486  { &kFitDmSq32, {hie * fabs(calc->GetDmsq32())} },
487  { thseed.var, thseed.seeds },
488  { &kFitDeltaInPiUnits,{0.,0.5, 1., 1.5} }
489  }, seedShifts);
490  minchi23= min(minchi23, thisminchi);
491 
492  if(corrSysts){
493  // Plot the systematic pulls and label with BFV
494  auto shifts = PlotSystShifts(auxShifts);
495  TString str = "Best fit " ;
496  for (auto &v:fitvars){
497  str += TString::Format(" %s=%.3f ",v->LatexName().c_str(),v->GetValue(calc));
498  }
499  str+= TString::Format(" LL=%.3f", thisminchi);
500  shifts->SetTitle(str);
501  gPad->Print(outdir + "debug_slice_shifts_bestfit_" + suffix +
502  (hie>0? "NH": "IH") + thseed.label + ".pdf");
503  }
504 
505  ResetOscCalcToDefault(calc);
506  auxShifts.ResetToNominal();
507  }
508  }
509  std::cout << "\nFound overall minchi " << minchi23 << std::endl;
510 
511  //We'll save to this file, make pretty plots later
512  TFile * outfile = new TFile(outfilename+".root","recreate");
513  TFile * outFCfile = new TFile(outFCfilename+".root","recreate");
514 
515  outfile->cd();
516  TVectorD v(1);
517  v[0] = minchi23;
518  v.Write("overall_minchi");
519 
520  //HACK HACK
521  int steps = 80;// 5;
522 
523  //Now create all the surfaces
524  for(int hie: {-1, +1}){
525 
526  ResetOscCalcToDefault(calc);
527  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
528 
529  std::cout << "Starting profile " << (hie>0? "NH ": "IH") << "\n\n";
530 
531  struct ProfDef{
532  const IFitVar * fitvar;
533  double minX;
534  double maxX;
535  std::vector <const IFitVar * > profvars;
536  std::vector <TString> profvarnames;
537  std::map <const IFitVar *, std::vector <double> >profseeds;
538  TString shortname;
539  };
540 
541  for(auto const & thseed:th23seeds){
542  ProfDef profdef;
543  const IFitVar* kCurTheta23 = thseed.var;
544 
545  if(!th23slice && !dmsqSlice){
546  profdef.fitvar = &kFitDeltaInPiUnits;
547  profdef.minX = 0;
548  profdef.maxX = 2;
549  profdef.profvars = {&kFitDmSq32, &kFitSinSq2Theta13, kCurTheta23};
550  profdef.profvarnames = {"DmSq32","SinSq2Theta13","SinSqTheta23"};
551  profdef.profseeds = {{kCurTheta23, thseed.seeds}};
552  profdef.shortname="delta";
553  }
554 
555  if(th23slice){
556  profdef.fitvar = &kFitSinSqTheta23;
557  profdef.minX = (nueOnly ? 0:0.3);
558  profdef.maxX = (nueOnly ? 1:0.7);
559  profdef.profvars = {&kFitDmSq32, &kFitSinSq2Theta13, &kFitDeltaInPiUnits};
560  profdef.profvarnames = {"DmSq32","SinSq2Theta13","DeltaCPpi"};
561  profdef.profseeds = {{&kFitDeltaInPiUnits,{0.5,1.,1.5}}};
562  profdef.shortname="th23";
563  }
564 
565  if(dmsqSlice){
566  profdef.fitvar = &kFitDmSq32;
567  profdef.minX = hie > 0 ? 2e-3 : -3e-3;
568  profdef.maxX = hie > 0 ? 3e-3 : -2e-3;
569  profdef.profvars = {kCurTheta23, &kFitSinSq2Theta13, &kFitDeltaInPiUnits};
570  profdef.profvarnames = {"SinSqTheta23","SinSq2Theta13","DeltaCPpi"};
571  profdef.profseeds = { {kCurTheta23, thseed.seeds},
572  {&kFitDeltaInPiUnits,{0.,0.5,1.,1.5}}, };
573  profdef.shortname="dmsq";
574  }
575 
576  std::map<const IFitVar*, TGraph*> profVarsMap;
577  std::map<const ISyst*, TGraph*> systMap;
578  MakeMaps (profdef.profvars, profVarsMap, systs, systMap);
579 
580  std::cout << " Profile "<< thseed.label << "\n\n";
581 
582  auto slice = SqrtProfile(exptThis, calc,
583  profdef.fitvar, steps, profdef.minX, profdef.maxX,
584  minchi23,
585  profdef.profvars,
586  systs,
587  profdef.profseeds,
588  seedShifts,
589  profVarsMap, systMap);
590  //Save!
591  TString hieroctstr = (hie > 0 ? "NH" : "IH") + thseed.label;
592 
593  outfile->cd();
594  slice->Write((TString ) "slice_" + profdef.shortname + "_" + hieroctstr);
595  SaveMaps (outFCfile, hieroctstr,
596  profdef.profvars, profdef.profvarnames, profVarsMap, systs, systMap);
597 
598  }//end seeds
599  }//end hie
600 
601  delete outfile;
602  std::cout << "\nProfiles saved to " << outfilename << std::endl;
603  return;
604  }//end createfile
605 
606 
607 
608 
609  //////////////////////////////////////////////////////////////////////
610  // Make plots
611  //////////////////////////////////////////////////////////////////////
612  else{
613  std::cout << "Making plots from " << outfilename+".root" << std::endl;
614  TFile * infile = new TFile (outfilename+".root","read");
615 
616  auto mins =* (TVectorD*)infile->Get("overall_minchi");
617  double minchi23 = mins[0];
618 
619  std::vector <TString> sliceNames = {};
620  if (!th23slice && ! dmsqSlice) sliceNames = {"slice_delta_"};
621  if (th23slice) sliceNames = {"slice_th23_"};
622  if (dmsqSlice) sliceNames = {"slice_dmsq_"};
623 
624  for(TString sliceName:sliceNames){
625  std::vector < std::pair <TGraph*, TString> > graphs;
626 
627  if(sliceName.Contains("delta")) {
628  for (TString hie:{"NH","IH"}){
629  for(auto const & thseed:th23seeds){
630  auto h = (TH1D*) infile ->Get( sliceName + hie + thseed.label);
631  if(!h) {
632  std::cerr << "Problem " << sliceName + hie + thseed.label << "\n";
633  exit(1);
634  }
635  graphs.push_back({HistoToCyclicGraph(h),hie + " " + thseed.label});
636  }
637  }
638  DrawSliceCanvas(sliceName + (fccorr?"fccorr":""), 5, 0, 2);
639  }
640  else {
641  for (TString hie:{"NH","IH"}) {
642  for(auto const & thseed:th23seeds){
643  auto h = (TH1D*) infile ->Get( sliceName + hie + thseed.label);
644  if(!h) {std::cerr << "Problem " << sliceName + hie + thseed.label << "\n"; }
645  graphs.push_back({new TGraph(h),hie + " " + thseed.label});
646  }
647  }
648 
649  if(sliceName.Contains("th23")) {
650  DrawSliceCanvas(sliceName + (fccorr?"fccorr":""), 3,
651  (nueOnly ? 0 : 0.32), (nueOnly ? 1 : 0.72));
652  }
653 
654  if(sliceName.Contains("dmsq")) {
655  DrawSliceCanvas(sliceName + (fccorr?"fccorr":""), 3, 2.15, 2.75);
656  }
657  }
658 
659  if(fccorr) {
660  if(!joint || !corrSysts || decomp != "combo")
661  { std::cerr << "\nWe don't have FC\n\n"; exit(1);}
662 
663  TString fcFolder="/nova/ana/nu_e_ana/Ana2017/FC/";
664  for (auto &gr:graphs) {
665  TString fcSliceN = "FCCol_";
666  fcSliceN += gr.second ;
667  if(th23slice) fcSliceN += "_slice_ssth23";
668  else if(dmsqSlice) fcSliceN += "_slice_dmsq32";
669  else fcSliceN += "_slice_delta";
670  fcSliceN.ReplaceAll (" LO","_LO"); //gah
671  fcSliceN.ReplaceAll (" UO","_UO");
672  fcSliceN.ReplaceAll (" ","");
673  bool useFourier = (fcSliceN.Contains("delta") &&
674  (fcSliceN.Contains("IH") /*|| fcSliceN.Contains("LO")*/)) ;
675  gr.first = FCCorrectSlice(gr.first, fcFolder + fcSliceN + ".root",
676  fourierFit, !useFourier, true);
677  }
678 
679  }//end fccorr
680 
681  for (auto &gr:graphs) {
682  if(gr.second.Contains("IH")) gr.first->SetLineColor(kPrimColorIH);
683  if(gr.second.Contains("NH")) gr.first->SetLineColor(k3SigmaConfidenceColorNH);
684  if(gr.second.Contains("UO")) gr.first->SetLineStyle(7);
685  gr.first->SetLineWidth(3);
686  if(dmsqSlice){ //HACK
687  for (int i =0; i < gr.first->GetN(); ++i){
688  gr.first->SetPoint(i,1000*fabs(gr.first->GetX()[i]),gr.first->GetY()[i]);
689  }
690  }
691  gr.first->Draw("l same");
692  }
693 
694  gPad->RedrawAxis();
695 
696  auto leg = SliceLegend(graphs, !dmsqSlice && !th23slice);
697  leg->Draw();
698 
699  POTtag();
700  //gPad->SetGrid();
701  for(TString ext: {".pdf",".root"}) {
702  gPad->Print(outdir + "slice_" + suffix +
703  (fccorr?"_fccorr":"") +
704  (fourierFit?"_smooth":"") +
705  ext);
706  }//end ext
707 
708  bool annotate = 1;
709  if(annotate) {
710 
711  for (auto &gr:graphs) {
712 
713  HighlightSignPoints(gr.first, sliceName);
714  }
715  gPad->SetGrid();
716  gPad->Print(outdir + "points_slice_" + suffix +
717  (fccorr?"_fccorr":"") +
718  (fourierFit?"_smooth":"") +
719  ".pdf");
720  }
721  }//end slicename
722 
723 
724  bool debug = true;
725 
726  if(debug){
727 
728  TFile * inFCfile = new TFile(outFCfilename+".root","read");
729  std::cout << "\n\nOpen profiled " << outFCfilename+".root\n" << std::endl;
730 
731  std::vector <TString> strprof = {"DmSq32","SinSq2Theta13","SinSqTheta23"};
732  if(th23slice) strprof = {"DmSq32","SinSq2Theta13","DeltaCPpi"};
733  if(dmsqSlice) strprof = {"SinSqTheta23","SinSq2Theta13","DeltaCPpi"};
734 
735  auto systs = GetJointFitSystematicList(corrSysts, !numuOnly, !nueOnly, true);
736 
737  for (const ISyst* s : systs){
738  strprof.push_back(s->ShortName());
739  }
740 
741  for(auto const & str:strprof){
742  new TCanvas;
743  double miny=-1.5, maxy=1.5;
744  if(str=="DeltaCPpi") {miny=0; maxy=2;}
745  if(str=="DmSq32") {miny=2.2; maxy=2.8;}
746  if(str=="SinSq2Theta13") {miny=0.080; maxy=0.086;}
747  if(str=="SinSqTheta23") {miny=0.4; maxy=0.6;}
748 
749  auto axes = new TH2F("",";#delta_{CP}/#pi;" + str, 100, 0, 2, 100, miny, maxy);
750  if(dmsqSlice) axes = new TH2F("",";|#Deltam^{2}_{32}| (10^{-3} eV^{2});" + str, 100, 2.2, 2.8, 100, miny, maxy);
751  if(th23slice) axes = new TH2F("",";sin^{2}#theta_{23};" + str, 100, 0.3, 0.7, 100, miny, maxy);
752  axes->Draw("hist");
753 
754  for (TString hie:{"NH","IH"}){
755  for(auto const & thseed:th23seeds){
756  auto gr = (TGraph*) inFCfile->Get(hie + thseed.label + "_" + str);
757  if(!gr) {
758  std::cerr << "Problem " << hie + thseed.label + "_" + str << "\n";
759  exit(1);
760  }
761  if(str=="DmSq32"){
762  for (int i = 0; i < (int) gr->GetN(); ++i){
763  gr->SetPoint(i, gr->GetX()[i], 1000* fabs(gr->GetY()[i]));
764  }
765  }
766  if(dmsqSlice){
767  for (int i = 0; i < (int) gr->GetN(); ++i){
768  gr->SetPoint(i, 1000*fabs(gr->GetX()[i]), gr->GetY()[i]);
769  }
770  }
771  gr->SetLineColor(hie=="NH"?k3SigmaConfidenceColorNH:kPrimColorIH);
772  if(thseed.label.Contains("UO")) gr->SetLineStyle(7);
773  gr->SetLineWidth(3);
774  gr->SetMarkerColor(gr->GetLineColor());
775  gr->Draw("l same");
776  }//end th23
777  }//end hie
778  gPad->SetGrid();
779  gPad->Print(outdir + "debug_slice_prof_" + suffix + "_"+ str + ".pdf");
780  }//end profile vars
781  }//end debug
782 
783  }//end plots
784 }//end all
Spectrum * GetFakeData(const IPrediction *pred, osc::IOscCalc *calc, const double pot, TH1D *cosmics=0)
const Color_t kPrimColorIH
Definition: Style.h:64
std::vector< std::pair< TH1D *, double > > GetNumuCosmics2017(const int nq=4, bool GetFromUPS=false)
void DrawSliceCanvas(TString slicename, const double ymax, const double xmin=0, const double xmax=2.)
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::map< std::string, double > xmax
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
TH1 * PlotSystShifts(const SystShifts &shifts, bool sortName)
Definition: Plots.cxx:1495
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
double delta
Definition: runWimpSim.h:98
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 MakeMaps(std::vector< const IFitVar * > profvars, std::map< const IFitVar *, TGraph * > &profMap, std::vector< const ISyst * > systs, std::map< const ISyst *, TGraph * > &systMap)
void XAxisDeltaCPLabels(TH1 *axes, bool t2kunits)
Label the x-axis with fractions of pi.
T sqrt(T number)
Definition: d0nt_math.hpp:156
const FitSinSqTheta23UpperOctant kFitSinSqTheta23UpperOctant
Definition: FitVars.cxx:16
void SaveMaps(TDirectory *out, std::map< std::string, IDecomp * > decomps_nominal, std::map< std::string, std::map< std::string, std::map< int, IDecomp * > > > decomps_shifted, std::map< std::string, PredictionNoExtrap * > predsNE_nominal, std::map< std::string, std::map< std::string, std::map< int, PredictionNoExtrap * > > > predsNE_shifted, std::map< std::string, PredictionSterile * > predsSt_nominal, std::map< std::string, std::map< std::string, std::map< int, PredictionSterile * > > > predsSt_shifted)
Save all of the objects in the input maps to the out directory/file.
Definition: PPFXHelper.h:1077
OStream cerr
Definition: OStream.cxx:7
const IPrediction * GetNuePrediction2017(std::string decomp, osc::IOscCalc *calc, bool corrSysts, bool GetFromUPS=false)
void abs(TH1 *hist)
const double kAna2017POT
Definition: Exposures.h:174
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
Float_t tmp
Definition: plot.C:36
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1481
static SystShifts Nominal()
Definition: SystShifts.h:34
osc::OscCalcDumb calc
TGraph * signFuncGr
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
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
virtual void SetDmsq32(const T &dmsq32)=0
void joint_fit_2017_slices(std::string decomp="prop", bool createFile=false, bool corrSysts=false, TString options="joint", bool th23slice=false, bool octantSlice=true, bool dmsqSlice=false, bool fccorr=false, bool fourierFit=false)
double SignificanceForUpValue(double up) const
Definition: FCBin.cxx:59
TGraph * FCCorrectSlice(TGraph *sqrtslice, TString fcFileName, bool fourierFit, bool useRoot=false, bool fccol=false)
const char * label
const XML_Char const XML_Char * data
Definition: expat.h:268
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.
Double_t ymax
Definition: plot.C:25
const XML_Char * s
Definition: expat.h:262
string outfilename
knobs that need extra care
const int nbins
Definition: cellShifts.C:15
std::pair< TH1D *, double > GetNueCosmics2017(bool GetFromUPS=false)
void SmoothWithFourierFit(TH1 *hist, int nbins=80, double minX=0, double maxX=2, int fourierFitN=3)
std::vector< const IPrediction * > GetNumuPredictions2017(const int nq=4, bool useSysts=true, bool GetFromUPS=false)
string infile
const Color_t k3SigmaConfidenceColorNH
Definition: Style.h:78
std::vector< Spectrum * > GetNumuData2017(const int nq=4)
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()
void DrawSignMarker(double x, double y, int color)
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
Float_t f1
virtual T GetDmsq32() const
Definition: IOscCalc.h:86
TLegend * SliceLegend(std::vector< std::pair< TGraph *, TString > > &graphs, bool isDelta)
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)
const XML_Char * prefix
Definition: expat.h:380
double sigma(TH1F *hist, double percentile)
TF1 * Fit() const
Definition: Utilities.cxx:694
void FindSignPoint(TGraph *gr, double sigma, double lowx=1, double highx=2, int color=kMagenta)
OStream cout
Definition: OStream.cxx:6
TGraph * HistoToCyclicGraph(TH1 *h)
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
Combine multiple component experiments.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const double kAna2017FullDetEquivPOT
Definition: Exposures.h:195
TDirectory * dir
Definition: macro.C:5
const NOvARwgtSyst kMECq0ShapeSyst2017("MECq0Shape","MEC q_{0} shape", novarwgt::kMECq0ShapeSyst2017)
Definition: MECSysts.h:41
const FitSinSqTheta23LowerOctant kFitSinSqTheta23LowerOctant
Definition: FitVars.cxx:17
void ResetToNominal()
Definition: SystShifts.cxx:141
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
exit(0)
void HighlightSignPoints(TGraph *gr, TString grname)
void POTtag()
Double_t signFunc(Double_t *x, Double_t *t)
And supporting functions:
assert(nhit_max >=nhit_nbins)
Interface definition for fittable variables.
Definition: IFitVar.h:16
c1
Definition: demo5.py:24
std::vector< std::pair< const ISyst *, const ISyst * > > GetCorrelations(bool isNue, bool isFHC, bool ptExtrap)
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
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
void AddCyclicPoints(TGraph *g)
Float_t e
Definition: plot.C:35
const Float_t kBlessedTitleFontSize
Definition: Style.h:89
#define for
Definition: msvc_pragmas.h:3
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)
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
static constexpr Double_t sr
Definition: Munits.h:164
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.
TGraph * SqrtProfile(const IExperiment *expt, osc::IOscCalcAdjustable *calc, const IFitVar *v, int nbinsx, double minx, double maxx, double minchi, std::vector< const IFitVar * > profVars, std::vector< const ISyst * > profSysts, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts, std::map< const IFitVar *, TGraph * > &profVarsMap, std::map< const ISyst *, TGraph * > &systsMap, MinuitFitter::FitOpts opts)
Forward to Profile but sqrt the result for a crude significance.
Definition: Fit.cxx:143
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17