joint_fit_2019_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 
21 //#include "./joint_fit_2019_loader_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 
39 #include <algorithm>
40 #include <vector>
41 #include <string>
42 
43 using namespace ana;
44 
45 
46 TGraph * signFuncGr;
47 Double_t signFunc(Double_t *x, Double_t *t);
48 void DrawSignMarker(double x, double y, int color);
49 void FindSignPoint(TGraph* gr, double sigma, double lowx, double highx, int color);
50 void HighlightSignPoints(TGraph* gr, TString grname);
51 void MakeMaps ( std::vector <const IFitVar*> profvars,
52  std::map<const IFitVar*, TGraph*> &profMap,
53  std::vector <const ISyst* > systs,
54  std::map<const ISyst*, TGraph*> &systMap);
55 void SaveMaps ( TDirectory * dir,
56  TString prefix,
57  std::vector <const IFitVar*> profvars,
58  std::vector <TString > profnames,
59  std::map<const IFitVar*, TGraph*> &profMap,
60  std::vector <const ISyst* > systs,
61  std::map<const ISyst*, TGraph*> &systMap);
62 void POTtag(bool both, bool FHCOnly, bool RHCOnly);
63 void AverageCyclicPoints(TGraph* g);
64 void AddCyclicPoints(TGraph* g);
65 TGraph* HistoToCyclicGraph(TH1* h, bool useRoot);
66 void SmoothWithFourierFit(TH1* hist,
67  int nbins, double minX, double maxX,
68  int fourierFitN);
69 //why is this function so ugly
70 TGraph * SmoothWithFourierFit(TGraph* gr0, bool isDelta,
71  bool useRoot, int fourierFitN);
72 TGraph* FCCorrectSlice(TGraph* sqrtslice, TString fcFileName,
73  bool fourierFit, std::string plot_name,
74  bool useRoot, bool fccol);
75 
76 
77 
78 void joint_fit_2019_slices(bool createFile = false,
79  bool corrSysts = true,
80  bool runOnGrid = false,
81  TString options="joint_realData_both",
82  bool th23slice = false,
83  bool octantSlice = true,
84  bool dmsqSlice = false,
85  bool th13Slice = false,
86  bool fccorr = false,
87  bool fourierFit = true)
88 
89 {
90  bool nueOnly = options.Contains("nueOnly");
91  bool numuOnly = options.Contains("numuOnly");
92  bool joint = options.Contains("joint");
93  assert (nueOnly || numuOnly || joint);
94 
95  bool FHCOnly = options.Contains("FHCOnly");
96  bool RHCOnly = options.Contains("RHCOnly");
97  bool both = options.Contains("both");
98  assert (FHCOnly || RHCOnly || both);
99 
100  bool fake2018 = options.Contains("fake2018");
101  bool realData = options.Contains("realData");
102 
103  auto suffix = options;
104  if(corrSysts) suffix+="_systs";
105  else suffix+="_stats";
106 
107  TString tag;
108  if(FHCOnly) tag="FHCOnly/";
109  if(RHCOnly) tag="RHCOnly/";
110  if(both) tag = "RHC_and_FHC/";
111 
112  if(th23slice) suffix+="_th23";
113  if(dmsqSlice) suffix+="_dmsq";
114  if(th13Slice) suffix+="_th13";
115  if(!octantSlice)suffix+="_noOct";
116 
117  // official path commented so things aren't being overwritten
118  TString outdir = "/nova/ana/nu_e_ana/Ana2019/Results/"+tag+"slices/FCInputs/";
119  if (runOnGrid) outdir = "./";
120 
121  TString outFCfilename (outdir + "slices_FCInput_2019_" + suffix);
122  if(joint && corrSysts && !th23slice && ! dmsqSlice && !th13Slice) outFCfilename = outdir + "slices_FCInput_2019_" +suffix;
123 
124  // outdir = "/nova/ana/nu_e_ana/Ana2019/Results/"+tag+"slices/";
125  outdir="./results/blessed_plots/slices/";
126  outdir += corrSysts?"syst/":"stat/";
127  if (runOnGrid) outdir = "./";
128 
129  TString outfilename (outdir + "hist_slices_2019_" + suffix);
130 
131  if(createFile){
132 
133  //////////////////////////////////////////////////
134  // Load Nue and Numu experiments
135  //////////////////////////////////////////////////
136  //need numu only for prestage seeds
137 
138  struct predictions {
139  const string name;
140  const IPrediction * pred;
141  std::pair <TH1D*, double> cos;
142  double pot;
143  };
144 
145  std::vector <predictions> preds;
146  std::vector <predictions> preds_numu_only;
147  std::vector <Spectrum * > data;
148  std::vector <Spectrum * > data_numu_only;
149  std::vector <const IExperiment * > expts;
150  std::vector <const IExperiment * > expts_numu_only;
151 
152  auto calc_fake = DefaultOscCalc();
153  if(fake2018) SetFakeCalc(calc_fake, 0.58, 2.51e-3, 0.17);
154  else if(!realData) {std::cerr << "need setting for data\n"; exit(1);}
155 
156  if(!numuOnly) {
157  if(FHCOnly || both ) {
158  preds.push_back({"nue_fhc", GetNuePrediction2019("combo", DefaultOscCalc(), corrSysts, "fhc", false, runOnGrid), GetNueCosmics2019("fhc", runOnGrid), GetPOT()});
159  }
160  if(RHCOnly || both ) {
161  preds.push_back({"nue_fhc", GetNuePrediction2019("prop", DefaultOscCalc(), corrSysts, "rhc", false, runOnGrid), GetNueCosmics2019("rhc", runOnGrid), GetPOT(false)});
162  }
163  }
164 
165  int nnumu = 4;
166  if(!nueOnly) {
167  if(FHCOnly || both ) {
168  auto numu_preds = GetNumuPredictions2019(nnumu, corrSysts, "fhc", runOnGrid);
169  auto numu_cosmics = GetNumuCosmics2019(nnumu, "fhc", runOnGrid);
170  for (int i = 0; i< nnumu; i++ ){
171  preds.push_back({"numu_fhc_"+std::to_string(i+1), numu_preds[i], numu_cosmics[i], GetPOT()});
172  preds_numu_only.push_back({"numu_fhc_"+std::to_string(i+1), numu_preds[i], numu_cosmics[i], GetPOT()});
173  }
174  }
175  if(RHCOnly || both ) {
176  auto numu_preds = GetNumuPredictions2019(nnumu, corrSysts, "rhc", runOnGrid);
177  auto numu_cosmics = GetNumuCosmics2019(nnumu, "rhc", runOnGrid);
178  for (int i = 0; i< nnumu; i++ ){
179  preds.push_back({"numu_rhc_"+std::to_string(i+1), numu_preds[i], numu_cosmics[i], GetPOT(false)});
180  preds_numu_only.push_back({"numu_rhc_"+std::to_string(i+1), numu_preds[i], numu_cosmics[i], GetPOT(false)});
181  }
182  }
183  }
184 
185 
186  if(realData){
187  if(!numuOnly){
188  if(FHCOnly || both ) data.push_back(GetNueData2019("fhc", runOnGrid));
189  if(RHCOnly || both ) data.push_back(GetNueData2019("rhc", runOnGrid));
190  }
191  if(!nueOnly){
192  if(FHCOnly || both ){
193  auto numu_data = GetNumuData2019(nnumu, "fhc", runOnGrid);
194  data.insert(data.end(),numu_data.begin(), numu_data.end());
195  data_numu_only.insert(data_numu_only.end(),numu_data.begin(), numu_data.end());}
196  if(RHCOnly || both ){
197  auto numu_data = GetNumuData2019(nnumu, "rhc", runOnGrid);
198  data.insert(data.end(),numu_data.begin(), numu_data.end());
199  data_numu_only.insert(data_numu_only.end(),numu_data.begin(), numu_data.end());}
200  }
201  }
202 
203  for(int i = 0; i < int(preds.size()); ++i){
204  double POT = preds[i].pot;
205  auto thisdata = GetFakeData (preds[i].pred, calc_fake, POT, preds[i].cos.first);
206  if(realData) thisdata = data[i];
207  cout<<i<<" "<<preds[i].name<<" POT "<<POT<<" tot MC "<<preds[i].pred->Predict(calc_fake).Integral(POT)<<" cos "<<preds[i].cos.first->Integral()<<" cos er "<<preds[i].cos.second<<" analyze data "<<thisdata->Integral(POT)<<endl;
208  expts.push_back(new SingleSampleExperiment(preds[i].pred, *thisdata, preds[i].cos.first, preds[i].cos.second));
209  }
210 
211  cout<<"Make numu only experiment to get the seeds"<<endl;
212  //make numu only experiment for seeds:
213 
214  for(int i = 0; i < (int) preds_numu_only.size(); ++i){
215  double POT = preds_numu_only[i].pot;
216  auto thisdata = GetFakeData (preds_numu_only[i].pred, calc_fake, POT, preds_numu_only[i].cos.first);
217  if(realData) thisdata = data_numu_only[i];
218  cout<<i<<" "<<preds_numu_only[i].name<<" POT "<<POT<<" tot MC "<<preds_numu_only[i].pred->Predict(calc_fake).Integral(POT)<<" cos "<<preds_numu_only[i].cos.first->Integral()<<" cos er "<<preds_numu_only[i].cos.second<<" analyze data "<<thisdata->Integral(POT)<<endl;
219  expts_numu_only.push_back(new SingleSampleExperiment(preds_numu_only[i].pred, *thisdata, preds_numu_only[i].cos.first, preds_numu_only[i].cos.second));
220  }
221  ////////////////////////////////////////////////////////////
222  // Add constraints, make experiments
223  ////////////////////////////////////////////////////////////
224  std::cout << "\nCreating multiexperiment\n" << std::endl;
225 
226  if (!th13Slice) {expts.push_back(WorldReactorConstraint2017()); std::cout << "Adding WorldReactorConstraint2017\n";}
227  else {std::cout << "No reactor constraint, that's th13 slice\n";}
228  if(nueOnly) {
229  std::cout << "Adding Dmsq32ConstraintPDG2017\n";
230  expts.push_back(&kDmsq32ConstraintPDG2017);
231  }
232  std::cout << "Creating Multiexperiment with a total of "
233  << expts.size() << " experiments\n\n" << std::endl;
234  auto exptThis = new MultiExperiment(expts);
235 
236  std::cout << "Creating Multiexperiment of numu only SimpleExp with a total of "
237  << expts_numu_only.size() << " experiments\n\n" << std::endl;
238  auto exptThis_numu_only = new MultiExperiment(expts_numu_only);
239 
240  ////////////////////////////////////////////////////////////
241  // Systematics
242  ////////////////////////////////////////////////////////////
243 
244  std::cout << "Systematics for the fit:\n";
245  auto systs = GetJointFitSystematicList(corrSysts, nueOnly, numuOnly, FHCOnly||both, RHCOnly||both, true);
246 
247  std::cout << "\n\nSystematic correlations...\n";
248  if(!nueOnly && ! numuOnly && corrSysts){
249  if(FHCOnly){
250  exptThis->SetSystCorrelations(0, GetCorrelations(true, true));
251  auto notfornumu = GetCorrelations(false, true);
252  for(int i =0; i < nnumu; ++i) exptThis->SetSystCorrelations(i+1, notfornumu);
253  }
254  if(RHCOnly){
255  exptThis->SetSystCorrelations(0, GetCorrelations(true, false));
256  auto notfornumu = GetCorrelations(false, true);
257  for(int i =0; i < nnumu; ++i) exptThis->SetSystCorrelations(i+1, notfornumu);
258  }
259  if(both){
260  exptThis->SetSystCorrelations(0, GetCorrelations(true, true));
261  exptThis->SetSystCorrelations(1, GetCorrelations(true, false));
262  auto notfornumu = GetCorrelations(false, true);
263  for(int i =0; i < 8; ++i) exptThis->SetSystCorrelations(i+2, notfornumu);
264  }
265  }
266  if (nueOnly && corrSysts){
267  if (FHCOnly) exptThis->SetSystCorrelations(0, GetCorrelations(true, true));
268  if (RHCOnly) exptThis->SetSystCorrelations(0, GetCorrelations(true, false));
269  if (both) {
270  exptThis->SetSystCorrelations(0, GetCorrelations(true, true));
271  exptThis->SetSystCorrelations(1, GetCorrelations(true, false));
272  }
273  }
274 
275  std::cout << "Systematics for the numu only fit:\n";
276  auto systs_numu_only = GetJointFitSystematicList(corrSysts, false, true, true, true, true);
277 
278  ////////////////////////////////////////////////////////////
279  // Fit
280  ////////////////////////////////////////////////////////////
281  std::cout << "Starting the fit" << std::endl;
282 
284 
285  SystShifts auxShifts = SystShifts::Nominal();
286 
287  // In case some systs need different seeds
288  std::vector <SystShifts> seedShifts = {};
289  if(corrSysts && !th23slice && !dmsqSlice && !th13Slice && octantSlice && (both || FHCOnly)){ //Cherenkov seed for deltaCP and neutron as well :(
290  cout<<"Add Cherenkov seeds \n";
291  for (double systshift:{-0.5, +0.5}){ //for the FHC Only or FHC+RHC ana
292  SystShifts tempShifts (&kAnaCherenkovSyst, systshift);
293  seedShifts.push_back(tempShifts);
294  }
295  if(both){
296  cout<<"Add neutron systematic seeds \n"; //For the joint RHC+FHC ana only
297  for (double systshift:{-0.5, +0.5}){
298  SystShifts tempShifts (&kNeutronVisEScalePrimariesSyst2018, systshift);
299  seedShifts.push_back(tempShifts);
300  }
301  }
302  }
303 
304 
305  //////////////////////////////////////////////////////////////////////
306  ///////////////////////// Seed controller ////////////////////////////
307  //////////////////////////////////////////////////////////////////////
308 
309  cout<<"------------------- Start prestage seeds --------------------------"<<endl;
310 
311  double minchi_numu = 1E20;
312  double pre_seed_th23;
313  double pre_seed_dmsq;
314  ResetOscCalcToDefault(calc);
315  auxShifts.ResetToNominal();
316 
317  double maxmix = 0.514; // from the numu results
318  double numu_pre_seedLONH, numu_pre_seedUONH, numu_pre_seedLOIH, numu_pre_seedUOIH, dmsq_numu_pre_seedNH, dmsq_numu_pre_seedIH;
319 
320  for(int hie:{1}){
321  std::cout << "\n\nFinding test best fit " << (hie>0? "NH " : "IH ") << "\n\n";
322  Fitter fitnumu_only(exptThis_numu_only, {&kFitSinSqTheta23, &kFitDmSq32}, {});
323 
324  auto thisminchi = fitnumu_only.Fit(calc, auxShifts ,
325  {{&kFitDmSq32,{hie*2.5e-3}},
326  {&kFitSinSqTheta23, {0.4} }}, {}, //seedShifts,
327  Fitter::kVerbose);
328 
329  if(thisminchi < minchi_numu) minchi_numu = thisminchi;
330 
331  pre_seed_th23 = kFitSinSqTheta23.GetValue(calc);
332  pre_seed_dmsq = kFitDmSq32.GetValue(calc);
333  }
334 
335  numu_pre_seedLONH = ((pre_seed_th23>maxmix)?(2*maxmix-pre_seed_th23):pre_seed_th23);
336  numu_pre_seedUONH = ((pre_seed_th23>maxmix)?pre_seed_th23:(2*maxmix-pre_seed_th23));
337 
338  ResetOscCalcToDefault(calc);
339  auxShifts.ResetToNominal();
340 
341  cout<<"------------------- End prestage seeds --------------------------"<<endl;
342 
343 
344  struct th23helper{
345  std::vector<double> seeds;
346  const IFitVar * var;
347  TString label;
348  };
349 
350 
351  std::vector <th23helper> th23seeds;
352 
353  th23helper temp;
354  temp.seeds = {0.3};
355  temp.var = kFitSinSqTheta23LowerOctant;
356  temp.label = "LO";
357 
358  //for NH:
359  if(octantSlice) {
360  th23seeds.push_back( { {0.499, numu_pre_seedLONH}, kFitSinSqTheta23LowerOctant, "LO"});
361  th23seeds.push_back( { {0.501, numu_pre_seedUONH}, kFitSinSqTheta23UpperOctant, "UO"});
362  }
363  else th23seeds.push_back({ {numu_pre_seedLONH, 0.5, numu_pre_seedUONH}, &kFitSinSqTheta23, ""});
364 
365  std::vector<double> delta_seeds = {0, 0.5, 1., 1.5};
366  std::vector<double> th23_old_seeds = {0.45, 0.5, 0.55};
367 
368  /////////////////////////////////////////////////////////////////////////
369  ////////////////////////// Different best fit searches //////////////////
370  /////////////////////////////////////////////////////////////////////////
371 
372  // Find the best fit points
373  double minchi23 = 1E20;
374  for(int hie:{-1, 1}){
375  for (auto & thseed:th23seeds){
376 
377  std::cout << "\n\nFinding best fit " << (hie > 0 ? "NH " : "IH ")
378  << thseed.label << ", "
379  << "ssth23 seeds ";
380  for (auto const & s:thseed.seeds) std::cout << s << ", ";
381  std::cout << std::endl;
382 
383  std::vector <const IFitVar*> fitvars = {&kFitDeltaInPiUnits,
384  thseed.var,
385  &kFitDmSq32,
387 
388  Fitter fit23(exptThis, fitvars, systs);
389  cout<<" pre dmsq seed is "<<pre_seed_dmsq<<endl;
390  auto thisminchi = fit23.Fit(calc, auxShifts,
391  {
392  { &kFitDmSq32, {hie*pre_seed_dmsq} },
393  { thseed.var, thseed.seeds },
394  { &kFitDeltaInPiUnits, delta_seeds }
395  }, seedShifts,
396  Fitter::kVerbose);
397  minchi23= min(minchi23, thisminchi);
398 
399  if(corrSysts){
400  // Plot the systematic pulls and label with BFV
401  auto shifts = PlotSystShifts(auxShifts);
402  TString str = "Best fit " ;
403  for (auto &v:fitvars){
404  str += TString::Format(" %s=%.3f ",v->LatexName().c_str(),v->GetValue(calc));
405  }
406  str+= TString::Format(" LL=%.6f", thisminchi);
407  shifts->SetTitle(str);
408  gPad->Update();
409  TLine *l=new TLine(gPad->GetUxmin(),0,gPad->GetUxmax(),0);
410  l->Draw("same");
411  gPad->Print(outdir + "debug_slice_shifts_bestfit_" + suffix +
412  (hie>0? "NH": "IH") + thseed.label + ".pdf");
413  }
414 
415  ResetOscCalcToDefault(calc);
416  auxShifts.ResetToNominal();
417  }
418  }
419  std::cout << "\nFound overall minchi " << minchi23 << "\n\n";
420 
421  cout<<"Check with oldstyle seeds "<<endl;
422  double minchi23test = 1E20;
423  for(int hie:{-1,1}){
424  std::cout << "\n\nFinding test best fit "
425  << (hie>0? "NH " : "IH ")
426  << "\n\n";
427  std::vector <const IFitVar*> fitvars = {&kFitDeltaInPiUnits,
429  &kFitDmSq32,
431 
432  Fitter fit23(exptThis, fitvars, systs);
433  auto thisminchi = fit23.Fit(calc,auxShifts ,{
434  {&kFitDmSq32,{hie*fabs(kFitDmSq32.GetValue(calc))}},
435  {&kFitSinSqTheta23, th23_old_seeds},
436  {&kFitDeltaInPiUnits, delta_seeds}
437  }, seedShifts
438  );
439  minchi23test= min(minchi23test, thisminchi);
440  if(corrSysts){
441  auto shifts = PlotSystShifts(auxShifts);
442  TString str = "Best fit " ;
443  for (auto &v:fitvars){
444  str += TString::Format(" %s=%.3f ",v->LatexName().c_str(),v->GetValue(calc));
445  }
446  str+= TString::Format(" LL=%.6f", thisminchi);
447  shifts->SetTitle(str);
448  gPad->Print(outdir + "debug_slice_shifts_bestfit_" + suffix +
449  (hie>0? "NH": "IH") + "test_chisq_with_no_octant_differentiation.pdf");
450  }
451 
452  ResetOscCalcToDefault(calc);
453  auxShifts.ResetToNominal();
454  }//end hie
455  std::cout << "\nFound overall test minchi " << minchi23test << std::endl;
456 
457 
458  cout<<"<<<<<<<<<<<<<<<<<<<<<<<<<< Test with only systematics free"<<endl;
459  minchi23test = 1E20;
460  ResetOscCalcToDefault(calc);
461  auxShifts.ResetToNominal();
462 
463  calc->SetdCP(0.16576*M_PI);
464  calc->SetTh23(asin(sqrt(0.584509)));
465  calc->SetDmsq32(2.50515e-3);
466  calc->SetTh13(asin(sqrt(0.0824663))/2);
467 
468  Fitter fit23(exptThis, {}, systs);
469  auto thisminchi = fit23.Fit(calc,auxShifts ,{}, seedShifts,
470  Fitter::kVerbose);
471  ofstream txtfile(outdir + "debug_slice_shifts_bestfit_" + suffix +
472  + "test_chisq_with_no_octant_test_with_osc_par_fixed.txt");
473  for (auto &v:systs){
474  txtfile<<v->ShortName().c_str()<<"="<<auxShifts.GetShift(v)<<endl;
475  }
476  txtfile.close();
477  minchi23test= min(minchi23test, thisminchi);
478  auto shifts = PlotSystShifts(auxShifts);
479  TString str = "Best fit " ;
480  str+= TString::Format(" LL=%.6f", thisminchi);
481  shifts->SetTitle(str);
482  gPad->Print(outdir + "debug_slice_shifts_bestfit_" + suffix +
483  + "test_chisq_with_no_octant_test_with_osc_par_fixed.pdf");
484  ResetOscCalcToDefault(calc);
485  auxShifts.ResetToNominal();
486  cout<<"<<<<<<<<<<<< This is not true BF : "<< minchi23test<<endl;
487 
488  /////////////////////////////////////////////////////////////////////////////
489  /////////////////////// Now, slices! ////////////////////////////////////////
490  ////////////////////////////////////////////////////////////////////////////
491 
492  //We'll save to this file, make pretty plots later
493  TFile * outfile = new TFile(outfilename+".root","recreate");
494  TFile * outFCfile = new TFile(outFCfilename+".root","recreate");
495 
496  outfile->cd();
497  TVectorD v(1);
498  v[0] = minchi23;
499  v.Write("overall_minchi");
500 
501  int steps = 60;// 80;
502 
503  //Now create all the slices
504  for(int hie: {-1, +1}){
505 
506  ResetOscCalcToDefault(calc);
507  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
508 
509  std::cout << "Starting profile " << (hie>0? "NH ": "IH") << "\n\n";
510 
511  struct ProfDef{
512  const IFitVar * fitvar;
513  double minX;
514  double maxX;
515  std::vector <const IFitVar * > profvars;
516  std::vector <TString> profvarnames;
517  std::map <const IFitVar *, std::vector <double> >profseeds;
518  TString shortname;
519  };
520 
521  for(auto const & thseed:th23seeds){
522  ProfDef profdef;
523  const IFitVar* kCurTheta23 = thseed.var;
524 
525  if(!th23slice && !dmsqSlice && !th13Slice){
526  profdef.fitvar = &kFitDeltaInPiUnits;
527  profdef.minX = 0;
528  profdef.maxX = 2;
529  profdef.profvars = {&kFitDmSq32, &kFitSinSq2Theta13, kCurTheta23};
530  profdef.profvarnames = {"DmSq32","SinSq2Theta13","SinSqTheta23"};
531  profdef.profseeds = {{kCurTheta23, thseed.seeds}, { &kFitDmSq32, {hie*pre_seed_dmsq}}, {&kFitSinSq2Theta13, {0.082}} };
532  profdef.shortname="delta";
533  }
534 
535 
536  if(th23slice){
537  profdef.fitvar = &kFitSinSqTheta23;
538  profdef.minX = (nueOnly ? 0:0.3);
539  profdef.maxX = (nueOnly ? 1:0.7);
540  profdef.profvars = {&kFitDmSq32, &kFitSinSq2Theta13, &kFitDeltaInPiUnits};
541  profdef.profvarnames = {"DmSq32","SinSq2Theta13","DeltaCPpi"};
542  profdef.profseeds = { {&kFitDeltaInPiUnits,delta_seeds}};
543  profdef.shortname="th23";
544  }
545 
546  if(th13Slice){
547  profdef.fitvar = &kFitSinSq2Theta13;
548  profdef.minX = 0.01;
549  profdef.maxX = 0.2;
550  profdef.profvars = {&kFitDmSq32, kCurTheta23, &kFitDeltaInPiUnits};
551  profdef.profvarnames = {"DmSq32","SinSqTheta23","DeltaCPpi"};
552  profdef.profseeds = { {kCurTheta23, thseed.seeds}, {&kFitDeltaInPiUnits,delta_seeds}, { &kFitDmSq32, {hie*pre_seed_dmsq}} };
553  profdef.shortname="th13";
554  }
555 
556  if(dmsqSlice){
557  profdef.fitvar = &kFitDmSq32;
558  profdef.minX = hie > 0 ? 2e-3 : -3e-3;
559  profdef.maxX = hie > 0 ? 3e-3 : -2e-3;
560  profdef.profvars = {kCurTheta23, &kFitSinSq2Theta13, &kFitDeltaInPiUnits};
561  profdef.profvarnames = {"SinSqTheta23","SinSq2Theta13","DeltaCPpi"};
562  profdef.profseeds = { {kCurTheta23, thseed.seeds}, {&kFitDeltaInPiUnits,delta_seeds}};
563  profdef.shortname="dmsq";
564  }
565 
566  std::map<const IFitVar*, TGraph*> profVarsMap;
567  std::map<const ISyst*, TGraph*> systMap;
568  MakeMaps (profdef.profvars, profVarsMap, systs, systMap);
569 
570  std::cout << " Profile "<< thseed.label << "\n";
571 
572  auto slice = SqrtProfile(exptThis, calc,
573  profdef.fitvar, steps, profdef.minX, profdef.maxX,
574  minchi23,
575  profdef.profvars,
576  systs,
577  profdef.profseeds,
578  seedShifts,
579  profVarsMap, systMap);
580  //Save!
581  TString hieroctstr = (hie > 0 ? "NH" : "IH") + thseed.label;
582 
583  outfile->cd();
584  slice->Write((TString ) "slice_" + profdef.shortname + "_" + hieroctstr);
585  SaveMaps (outFCfile, hieroctstr,
586  profdef.profvars, profdef.profvarnames, profVarsMap, systs, systMap);
587 
588  }//end seeds
589  }//end hie
590 
591  delete outfile;
592  std::cout << "\nProfiles saved to " << outfilename << std::endl;
593  return;
594  }//end createfile
595 
596  //////////////////////////////////////////////////////////////////////
597  // Make plots
598  //////////////////////////////////////////////////////////////////////
599  else{
600  TString infilename = "/nova/ana/nu_e_ana/Ana2019/Results/RHC_and_FHC/slices/";
601  infilename += corrSysts?"syst/":"stat/";
602  infilename += "hist_slices_2019_";
603  infilename += options;
604  infilename += corrSysts?"_systs":"_stats";
605  if(th23slice) infilename+="_th23";
606  if(dmsqSlice) infilename+="_dmsq";
607  if(th13Slice) infilename+="_th13";
608  if(!octantSlice)infilename+="_noOct";
609  infilename += ".root";
610 
611  std::cout << "Making plots from " << infilename << std::endl;
612  TFile * infile = new TFile (infilename,"read");
613 
614  auto mins =* (TVectorD*)infile->Get("overall_minchi");
615  double minchi23 = mins[0];
616 
617  std::vector <TString> sliceNames = {};
618  if (!th23slice && ! dmsqSlice && !th13Slice) sliceNames = {"slice_delta_"};
619  if (th23slice) sliceNames = {"slice_th23_"};
620  if (dmsqSlice) sliceNames = {"slice_dmsq_"};
621  if (th13Slice) sliceNames = {"slice_th13_"};
622 
623  // it shouldn't play any role, but copy old style seeding here for sake of naming, everywhere is just .lable
624 
625  struct th23helper{
626  std::vector<double> seeds;
627  const IFitVar * var;
628  TString label;
629  };
630  std::vector <th23helper> th23seeds;
631  if(octantSlice) {
632  th23seeds.push_back( { {0.45}, kFitSinSqTheta23LowerOctant, "LO"});
633  th23seeds.push_back( { {0.55}, kFitSinSqTheta23UpperOctant, "UO"});
634  }
635  else th23seeds.push_back({ {0.45, 0.5, 0.55}, &kFitSinSqTheta23, ""});
636 
637  //end of old-style-seeding
638 
639  for(TString sliceName:sliceNames){
640  std::vector < std::pair <TGraph*, TString> > graphs;
641 
642  if(sliceName.Contains("delta")) {
643  for (TString hie:{"NH","IH"}){
644  for(auto const & thseed:th23seeds){
645  auto h = (TH1D*) infile ->Get( sliceName + hie + thseed.label);
646  if(!h) {
647  std::cerr << "Problem " << sliceName + hie + thseed.label << "\n";
648  exit(1);
649  }
650  if(fccorr)
651  graphs.push_back({new TGraph(h), hie + " " + thseed.label});
652  else
653  graphs.push_back({HistoToCyclicGraph(h, true),hie + " " + thseed.label});
654  }
655  }
656  DrawSliceCanvas(sliceName + (fccorr?"fccorr":""), 5, 0, 2);
657  }
658  else {
659  for (TString hie:{"NH","IH"}) {
660  for(auto const & thseed:th23seeds){
661  auto h = (TH1D*) infile ->Get( sliceName + hie + thseed.label);
662  if(!h) {std::cerr << "Problem " << sliceName + hie + thseed.label << "\n"; }
663  graphs.push_back({new TGraph(h),hie + " " + thseed.label});
664  }
665  }
666 
667  if(sliceName.Contains("th23")) {
668  DrawSliceCanvas(sliceName + (fccorr?"fccorr":""), 3,
669  (nueOnly ? 0 : 0.32), (nueOnly ? 1 : 0.72));
670  }
671 
672  if(sliceName.Contains("dmsq")) {
673  DrawSliceCanvas(sliceName + (fccorr?"fccorr":""), 3, 2.15, 2.75);
674  }
675  if(sliceName.Contains("th13")) {
676  DrawSliceCanvas(sliceName + (fccorr?"fccorr":""), 3, 0.05, 0.2);
677  }
678  }
679 
680  TString fcFolder="/nova/ana/nu_e_ana/Ana2019/FC/slices/";
681  std::string plot_name = "";
682  if(th23slice) {
683  plot_name = "ssth23";
684  fcFolder += "th23/";
685  }
686  else if(th13Slice) {
687  plot_name = "th13";
688  fcFolder += "th13/";
689  }
690  else if(dmsqSlice) {
691  plot_name = "dmsq";
692  fcFolder += "dmsq/";
693  }
694  else {
695  plot_name = "delta";
696  fcFolder += "delta/";
697  }
698 
699  if(fccorr) {
700  if(!joint || !corrSysts)
701  { std::cerr << "\nWe don't have FC\n\n"; exit(1);}
702 
703  // TString fcFolder="./FC/";
704  for (auto &gr:graphs) {
705  TString fcSliceN = plot_name;
706  fcSliceN += gr.second.Contains("NH")? "_nh" : "_ih";
707  if(!th23slice)
708  fcSliceN += gr.second.Contains("UO")? "uo" : "lo";
709  fcSliceN += ".root";
710 
711  // only use fourier fit for period boundary conditions
712  // Fourier fit doesn't like delta_UO_nh. Would need
713  // lots of harmonics to fit it correctly anyway
714  bool useFourier = (fcSliceN.Contains("delta") &&
715  !fcSliceN.Contains("delta_nhuo"));
716 
717  gr.first = FCCorrectSlice(gr.first, fcFolder + fcSliceN,
718  fourierFit, plot_name, !useFourier, true);
719  }
720 
721  }//end fccorr
722 
723  for (auto &gr:graphs) {
724  //if(gr.second.Contains("IH")) gr.first->SetLineColor(kPrimColorIH);
725  //if(gr.second.Contains("NH")) gr.first->SetLineColor(k3SigmaConfidenceColorNH);
726  if(gr.second.Contains("IH")) gr.first->SetLineColor(cih_dark);
727  if(gr.second.Contains("NH")) gr.first->SetLineColor(cnh);
728  if(gr.second.Contains("LO")) gr.first->SetLineStyle(7);
729  gr.first->SetLineWidth(3);
730  if(dmsqSlice){
731  for (int i =0; i < gr.first->GetN(); ++i){
732  gr.first->SetPoint(i,1000*fabs(gr.first->GetX()[i]),gr.first->GetY()[i]);
733  }
734  }
735  gr.first->Draw("l same");
736  }
737 
738  POTtag(both, FHCOnly, RHCOnly);
739  PreliminarySide();
740  gPad->RedrawAxis();
741 
742  auto leg = SliceLegend(graphs, !dmsqSlice && !th23slice && !th13Slice, 0.7, dmsqSlice);
743  leg->Draw();
744 
745  // outdir = "test_delta_slice/";
746  //gPad->SetGrid();
747  for(TString ext: {".pdf",".root"}) {
748  gPad->Print(outdir + "slice_" + suffix +
749  // gPad->Print(TString("./test/") + [>(corrSysts?"syst/":"stat/")+<] "slice_" + suffix +
750  (fccorr?"_fccorr":"") +
751  (fourierFit?"_smooth":"") +
752  ext);
753  }//end ext
754 
755  bool annotate = 1;
756  if(annotate) {
757 
758  for (auto &gr:graphs) {
759 
760  HighlightSignPoints(gr.first, sliceName);
761  }
762  gPad->SetGrid();
763  gPad->Print(outdir + "points_slice_" + suffix +
764  // gPad->Print(TString("./test/") + [>(corrSysts?"syst/":"stat/")+<] "points_slice_" + suffix +
765  (fccorr?"_fccorr":"") +
766  (fourierFit?"_smooth":"") +
767  ".pdf");
768  }
769  }//end slicename
770 
771 
772  bool debug = false;
773 
774  if(debug){
775 
776  TFile * inFCfile = new TFile(outFCfilename+".root","read");
777  std::cout << "\n\nOpen profiled " << outFCfilename+".root\n" << std::endl;
778 
779  std::vector <TString> strprof = {"DmSq32","SinSq2Theta13","SinSqTheta23"};
780  if(th23slice) strprof = {"DmSq32","SinSq2Theta13","DeltaCPpi"};
781  if(th13Slice) strprof = {"DmSq32","DeltaCPpi", "SinSqTheta23"};
782  if(dmsqSlice) strprof = {"SinSqTheta23","SinSq2Theta13","DeltaCPpi"};
783 
784  auto systs = GetJointFitSystematicList(corrSysts, numuOnly, nueOnly, true, true, true);
785 
786  for (const ISyst* s : systs){
787  strprof.push_back(s->ShortName());
788  }
789 
790  for(auto const & str:strprof){
791  new TCanvas;
792  double miny=-1.5, maxy=1.5;
793  if(str=="DeltaCPpi") {miny=0; maxy=2;}
794  if(str=="DmSq32") {miny=2.2; maxy=2.9;}
795  if(str=="SinSq2Theta13") {miny=0.080; maxy=0.086;}
796  if(str=="SinSqTheta23") {miny=0.3; maxy=0.8;}
797 
798  auto axes = new TH2F("",";#delta_{CP}/#pi;" + str, 100, 0, 2, 100, miny, maxy);
799  if(dmsqSlice) axes = new TH2F("",";|#Deltam^{2}_{32}| (10^{-3} eV^{2});" + str, 100, 2.2, 2.8, 100, miny, maxy);
800  if(th23slice) axes = new TH2F("",";sin^{2}#theta_{23};" + str, 100, 0.3, 0.7, 100, miny, maxy);
801  if(th13Slice) axes = new TH2F("",";sin^{2}#theta_{13};" + str, 100, 0.01, 0.2, 100, miny, maxy);
802  axes->Draw("hist");
803 
804  for (TString hie:{"NH","IH"}){
805  for(auto const & thseed:th23seeds){
806  auto gr = (TGraph*) inFCfile->Get(hie + thseed.label + "_" + str);
807  if(!gr) {
808  std::cerr << "Problem " << hie + thseed.label + "_" + str << "\n";
809  exit(1);
810  }
811  if(str=="DmSq32"){
812  for (int i = 0; i < (int) gr->GetN(); ++i){
813  gr->SetPoint(i, gr->GetX()[i], 1000* fabs(gr->GetY()[i]));
814  }
815  }
816  if(dmsqSlice){
817  for (int i = 0; i < (int) gr->GetN(); ++i){
818  gr->SetPoint(i, 1000*fabs(gr->GetX()[i]), gr->GetY()[i]);
819  }
820  }
821  //gr->SetLineColor(hie=="NH"?k2SigmaConfidenceColorNH:kPrimColorIH);
822  gr->SetLineColor(hie=="NH"?cnh:cih_dark);
823  if(thseed.label.Contains("LO")) gr->SetLineStyle(7);
824  gr->SetLineWidth(3);
825  gr->SetMarkerColor(gr->GetLineColor());
826  gr->Draw("l same");
827  }//end th23
828  }//end hie
829  gPad->SetGrid();
830  gPad->Print(outdir + "debug_slice_prof_" + suffix + "_"+ str + ".pdf");
831  }//end profile vars
832  }//end debug
833 
834  }//end plots
835 }//end all
836 
837 
838 /// And supporting functions:
839 
840 Double_t signFunc(Double_t *x, Double_t *t){
841 
842  double y = signFuncGr->Eval(x[0]);
843  if(t[0] > 0 && abs(y - t[0]) > 0.05) return 999.;
844  return abs(y - t[0]);
845 }
846 void DrawSignMarker(double x, double y, int color)
847 {
848  TMarker * m = new TMarker (x, y,29);
849  m->SetMarkerColor(color);
850  m->SetMarkerSize(3);
851  TString pstr = TString::Format("#bf{(%.2f,%.2f)}",x,y);
852  TLatex * ltx = new TLatex(x,y,pstr.Data());
853  ltx->SetTextAngle(45);
854  ltx->SetTextColor(color);
855  ltx->Draw();
856  m->Draw();
857 }
858 void FindSignPoint(TGraph* gr, double sigma, double lowx=1, double highx=2, int color = kMagenta)
859 {
860  double minX = gr->GetX()[0];
861  double maxX = gr->GetX()[gr->GetN()];
862  signFuncGr = gr;
863  TF1 * f1 = new TF1("sign",signFunc,minX,maxX,1);
864  f1->SetParameter(0,sigma);
865  double xmin = f1->GetMinimumX(lowx,highx);
866 
867  if(f1->Eval(xmin) < 5) DrawSignMarker(xmin, gr->Eval(xmin), color);
868 }
869 
870 
871 
872 void HighlightSignPoints(TGraph* gr, TString grname)
873 {
874  struct SPoint{
875  double sigma;
876  double xlow;
877  double xhigh;
878  };
879  std::vector <SPoint > sig_range;
880  //fish for local minima points / 123 sgma points
881  if(grname.Contains("dmsq")){
882  sig_range = {{0.,2.3,2.6},
883  {1,2.3,2.4},{1,2.4,2.6},
884  {2,2.2,2.5},{2,2.5,2.7}};
885  }
886  if(grname.Contains("th23")){
887  sig_range = {{0.,0,0.5},{0,0.5,0.7},{0,0.5,0.5001},
888  {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},
889  {2,0.3,0.5},{2,0.5,0.7}};
890  }
891  if(grname.Contains("th13")){
892  sig_range = {{0.,0,0.1},{0,0.1,0.15},{0,0.15,0.2},
893  {1,0.0,0.1},{1,0.1,0.15},{1,0.15,0.2},
894  {2,0.0,0.1},{2,0.1,0.2}};
895  }
896  if(grname.Contains("delta")){
897  sig_range = {{0.,1., 2.}, {0.,0., 1.},
898  {1., 0, 0.5},{1.,0.5, 1.},{1.,1, 1.7},{1.,1.7, 2.},
899  {2., 0, 0.5},{2.,0.5, 1.},{2.,1, 1.3},{2.,1.5, 2.},
900  {3., 0, 0.5}, {3,0.5,1.5}, {3,1.7,2.}};
901 
902  }
903  int color = kMagenta;
904  for (auto &sr:sig_range){
905  if(sr.sigma == 1 ) color = kGreen + 1;
906  if(sr.sigma == 2 ) color = kAzure ;
907  if(sr.sigma == 3 ) color = kOrange ;
908  FindSignPoint(gr,sr.sigma,sr.xlow, sr.xhigh, color);
909  }
910 }
911 
912 
913 void MakeMaps ( std::vector <const IFitVar*> profvars,
914  std::map<const IFitVar*, TGraph*> &profMap,
915  std::vector <const ISyst* > systs,
916  std::map<const ISyst*, TGraph*> &systMap)
917 {
918 
919  for (const IFitVar * var:profvars)
920  profMap.insert(std::pair<const IFitVar*, TGraph*> (var, new TGraph()));
921  for (const ISyst* syst : systs)
922  systMap.insert(std::pair<const ISyst*, TGraph*> (syst, new TGraph()));
923 }
924 
925 void SaveMaps ( TDirectory * dir,
926  TString prefix,
927  std::vector <const IFitVar*> profvars,
928  std::vector <TString > profnames,
929  std::map<const IFitVar*, TGraph*> &profMap,
930  std::vector <const ISyst* > systs,
931  std::map<const ISyst*, TGraph*> &systMap)
932 {
933  TDirectory *tmp = gDirectory;
934  dir->cd();
935  for (int i = 0; i < (int) profvars.size(); ++i){
936  profMap[profvars[i]]->Write((prefix + "_" + profnames[i]));
937  }
938  for (const ISyst* s : systs)
939  systMap[s]->Write((prefix + "_" + s->ShortName()));
940  tmp->cd();
941 }
942 
943 void POTtag(bool both, bool FHCOnly, bool RHCOnly)
944 {
945  TLatex* ltx = new TLatex();
946  auto ltx1 = new TLatex(.15, 0.93, TString::Format("NOvA FD"));
947  if (FHCOnly){ ltx = new TLatex(.42, 0.93, TString::Format("%g#times10^{20} POT equiv. of #nu data only", kAna2017FullDetEquivPOT/1E20).Data());}
948  if (RHCOnly){ ltx = new TLatex(.52, 0.93, TString::Format("%2.4g#times10^{20} POT of #bar{#nu} data only", kAna2019RHCPOT/1E20).Data());}
949  if(both) {ltx = new TLatex(.32, 0.93, TString::Format("%g#times10^{20} POT equiv #nu + %2.4g#times10^{20} POT #bar{#nu}", kAna2017FullDetEquivPOT/1E20, kAna2019RHCPOT/double(1E20)).Data());}
950  ltx1->SetNDC();
951  ltx1->SetTextAlign(11);
952  ltx1->SetTextSize(0.8*kBlessedLabelFontSize);
953  ltx1->Draw();
954  ltx->SetNDC();
955  ltx->SetTextAlign(11);
956  ltx->SetTextSize(0.8*kBlessedLabelFontSize);
957  ltx->Draw("same");
958 }
959 
960 void AverageCyclicPoints(TGraph* g)
961 {
962  double x0, y0;
963  double x2pi, y2pi;
964  g->GetPoint(g->GetN()-1, x2pi, y2pi);
965  g->GetPoint(1, x0, y0);
966  std::cout << "(x0, y0) = (" << x0 << ", " << y0 << ")\n";
967  std::cout << "(x2pi, y2pi) = (" << x2pi << ", " << y2pi << ")\n";
968 
969  double avgY = ( y0 + y2pi ) / 2;
970  std::cout << "Average = " << avgY << std::endl;
971  g->SetPoint(0, 0, avgY);
972  g->SetPoint(g->GetN(), 2, avgY);
973 
974  std::cout << "\nFound average cyclic point: (0, " << avgY << ")\n";
975  std::cout << "\nAdded average cyclic points " << 0 << ": " << g->GetX()[0] << " " << g->GetY()[0];
976  std::cout << ", " << g->GetN() << ": " << g->GetX()[g->GetN()-1] << " " << g->GetY()[g->GetN()-1] << "\n\n";
977 }
978 
979 
980 
981 void AddCyclicPoints(TGraph* g)
982 {
983  double x, y;
984  g->GetPoint(g->GetN()-1, x, y);
985  g->SetPoint(0, x-2, y);
986  g->GetPoint(1, x, y);
987  g->SetPoint(g->GetN(), x+2, y);
988 
989  std::cout << "\n Added cyclic points " << 0 << " " << g->GetX()[0]<< " " << g->GetY()[0]
990  << ", " << g->GetN()-1 << x+2 << y << "\n\n";
991 }
992 
993 
994 TGraph* HistoToCyclicGraph(TH1* h, bool useRoot = false)
995 {
996  TGraph* g = new TGraph;
997  g->SetPoint(0, 0, 0); // temporary
998  for(int i = 1; i <= h->GetNbinsX(); ++i)
999  g->SetPoint(i, h->GetXaxis()->GetBinCenter(i), h->GetBinContent(i));
1000 
1001  g->SetLineColor(h->GetLineColor());
1002  g->SetLineWidth(2);
1003 
1004  if(useRoot)
1006  else
1007  AddCyclicPoints(g);
1008 
1009  return g;
1010 }
1011 
1012 
1014  int nbins = 80, double minX = 0, double maxX = 2,
1015  int fourierFitN=3)
1016 {
1017  FitToFourier fitF (hist, minX, maxX, fourierFitN);
1018  auto f1= fitF.Fit();
1019  hist->Reset();
1020  hist->Add(f1);
1021 }
1022 //why is this function so ugly
1023 TGraph * SmoothWithFourierFit(TGraph* gr0,
1024  bool isDelta = true,
1025  bool useRoot = false,
1026  int fourierFitN = 2)
1027 {
1028  int nbins = gr0->GetN();
1029  double minX = gr0->GetX()[0];
1030  double maxX = gr0->GetX()[nbins-1];
1031  //TGraph to center bin
1032  minX = minX - (gr0->GetX()[1] - gr0->GetX()[0])/2;
1033  maxX = maxX + (gr0->GetX()[1] - gr0->GetX()[0])/2;
1034 
1035  //Fourier fit relies on periodicity, ignore extra points
1036  if(isDelta && !useRoot){
1037  nbins = 60;
1038  minX = 0;
1039  maxX = 2;
1040  }
1041  std::cout << "Smooth " << nbins << " nbins " << minX << ", " << maxX << "\n";
1042 
1043  TH1D * hist = new TH1D(UniqueName().c_str(),"",nbins,minX,maxX);
1044 
1045  for (int i = 0; i < nbins; ++i){
1046  double val = gr0->Eval(hist->GetBinCenter(i+1));
1047  hist->SetBinContent(i + 1, val * val); //use dchisq not chi2 to smooth
1048  }
1049  if(useRoot) hist->Smooth(2);
1050  else SmoothWithFourierFit(hist, nbins, minX, maxX, fourierFitN);
1051  //now return sqrt
1052  for (int i = 1; i <= nbins; ++i){
1053  hist->SetBinContent(i, sqrt(hist->GetBinContent(i)));
1054  }
1055  if(isDelta/* && !useRoot*/)
1056  return HistoToCyclicGraph(hist, useRoot);
1057  else return new TGraph(hist);
1058 }
1059 
1060 TGraph* FCCorrectSlice(TGraph* sqrtslice, TString fcFileName,
1061  bool fourierFit, std::string plot_name,
1062  bool useRoot = false, bool fccol=false)
1063 {
1064  bool th23plot = fcFileName.Contains("th23");
1065  bool dmsqplot = fcFileName.Contains("dmsq");
1066  bool dmsqnhuo = fcFileName.Contains("dmsq_nhuo");
1067  TString hiername = fcFileName.Contains("nh")? "NH":"IH";
1068 
1069  double minX = 0;
1070  double maxX = 2;
1071  if (th23plot) {minX = 0.3; maxX = 0.7;}
1072  if (dmsqplot) {
1073  minX = hiername=="NH" ? 2e-3 : -3e-3;
1074  maxX = hiername=="NH" ? 3e-3 : -2e-3;
1075  }
1076 
1077  //We want slice, not sqrt
1078  auto n= sqrtslice->GetN();
1079 
1080  auto sliceXH = new TGraph(n);
1081  for (int i = 0; i < n;i++ ){
1082  sliceXH->SetPoint(i,sqrtslice->GetX()[i],util::sqr(sqrtslice->GetY()[i]));
1083  }
1084  FCSurface *fcXH = 0;
1085  if(!fccol) fcXH = FCSurface::LoadFromFile(fcFileName.Data()).release();
1086  else {
1087  auto fccol = FCCollection::LoadFromFile(fcFileName.Data()).release();
1088  //Only one Y bin required for slices.
1089  // All slices this year had 60 bins. Need to find a better way
1090  // next year
1091  fcXH = new FCSurface(60,minX,maxX,1,0,1);
1092  fcXH->Add(*fccol, plot_name);
1093  }
1094 
1095 
1096  TGraph* sigXH = new TGraph;
1097  TGraph* pcXH = new TGraph;
1098 
1099  TFile* bestfitfile = new TFile("/nova/ana/nu_e_ana/Ana2019/Results/BestFits/bestfits_joint_realData_both_systs.root", "read");
1100  auto calcbf = (osc::IOscCalcAdjustable*)(LoadFrom<osc::IOscCalc>(bestfitfile, "NH_UO_calc")).release();
1101  double dmsq32bf = kFitDmSq32.GetValue(calcbf);
1102  for(int i = 0; i < n; ++i){
1103  const double delta = sliceXH->GetX()[i];
1104  const double upXH = sliceXH->GetY()[i];
1105 
1106  // Fall back to gaussian up if dchisq is too large or FCBin is empty (mostly for deltaCP IH slices)
1107  if(upXH > 13 || fcXH->GetBin(i+1,1)->Empty()){
1108  std::cout << " Falling back on gaussian up value: " << upXH
1109  << " for bin " << i << " , delta = " << delta << std::endl;
1110  sigXH->SetPoint(sigXH->GetN(), delta, sqrt(upXH));
1111  }
1112  else {
1113  double pXH;
1114  if(upXH > 0) pXH = fcXH->GetBin(i+1, 1)->SignificanceForUpValue(upXH);
1115  else pXH = 0; //HACK HACK
1116  sigXH->SetPoint(sigXH->GetN(), delta, PValueToSigma(1-pXH));
1117  }
1118  if(dmsqnhuo && (i+1 < n) && (delta < dmsq32bf) && (dmsq32bf < sliceXH->GetX()[i+1])){
1119  sigXH->SetPoint(sigXH->GetN(), dmsq32bf, 0.);
1120  }
1121  }
1122  if(fourierFit) {
1123  sigXH = SmoothWithFourierFit(sigXH, fcFileName.Contains("delta"), useRoot);
1124  }
1125 
1126  return sigXH;
1127 }
1128 
1129 // th23
1130 // cafe joint_fit_2019_slices.C false true false joint_realData_both true false false false true true
1131 
1132 // dmsq
1133 // cafe joint_fit_2019_slices.C false true false joint_realData_both false true true false true true
1134 
1135 // delta
1136 // cafe joint_fit_2019_slices.C false true false joint_realData_both false true false false true true
Spectrum * GetFakeData(const IPrediction *pred, osc::IOscCalc *calc, const double pot, TH1D *cosmics=0)
void DrawSliceCanvas(TString slicename, const double ymax, const double xmin=0, const double xmax=2.)
void joint_fit_2019_slices(bool createFile=false, bool corrSysts=true, 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)
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
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 GetValue(const osc::IOscCalcAdjustable *osc) const override
void PreliminarySide()
Definition: rootlogon.C:114
double delta
Definition: runWimpSim.h:98
TGraph * HistoToCyclicGraph(TH1 *h, bool useRoot)
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:18
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
std::vector< const IPrediction * > GetNumuPredictions2019(const int nq=4, bool useSysts=true, std::string beam="fhc", bool GetFromUPS=false, ENu2018ExtrapType numuExtrap=kNuMu, bool minimizeMemory=false, bool NERSC=false)
double maxy
const Dmsq32Constraint kDmsq32ConstraintPDG2017(2.45e-3, 0.05e-3, 2.52e-3, 0.05e-3)
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
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
void MakeMaps(std::vector< const IFitVar * > profvars, std::map< const IFitVar *, TGraph * > &profMap, std::vector< const ISyst * > systs, std::map< const ISyst *, TGraph * > &systMap)
virtual void SetTh13(const T &th13)=0
OStream cerr
Definition: OStream.cxx:7
void SmoothWithFourierFit(TH1 *hist, int nbins, double minX, double maxX, int fourierFitN)
void abs(TH1 *hist)
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
const NeutronVisEScaleSyst2018 kNeutronVisEScalePrimariesSyst2018(false)
void FindSignPoint(TGraph *gr, double sigma, double lowx, double highx, int color)
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
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 IPrediction * GetNuePrediction2019(std::string decomp, osc::IOscCalc *calc, bool corrSysts, std::string beam, bool isFakeData, bool GetFromUPS=false, bool minimizeMemory=false, bool NERSC=false)
virtual void SetDmsq32(const T &dmsq32)=0
#define M_PI
Definition: SbMath.h:34
double SignificanceForUpValue(double up) const
Definition: FCBin.cxx:59
TGraph * signFuncGr
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
TGraph * FCCorrectSlice(TGraph *sqrtslice, TString fcFileName, bool fourierFit, std::string plot_name, bool useRoot, bool fccol)
double PValueToSigma(double p)
Compute the equivalent number of gaussian sigma for a p-value.
void DrawSignMarker(double x, double y, int color)
const XML_Char * s
Definition: expat.h:262
string outfilename
knobs that need extra care
const int nbins
Definition: cellShifts.C:15
const DummyAnaSyst kAnaCherenkovSyst("Cherenkov","Cherenkov")
string infile
T GetShift(const ISyst *syst) const
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)
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 HighlightSignPoints(TGraph *gr, TString grname)
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
void AverageCyclicPoints(TGraph *g)
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
Spectrum * GetNueData2019(std::string beam, bool GetFromUPS=false)
const double kAna2019RHCPOT
Definition: Exposures.h:224
Combine multiple component experiments.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const double kAna2017FullDetEquivPOT
Definition: Exposures.h:195
std::vector< std::pair< TH1D *, double > > GetNumuCosmics2019(const int nq=4, std::string beam="fhc", bool GetFromUPS=false, bool NERSC=false)
void AddCyclicPoints(TGraph *g)
std::vector< Spectrum * > GetNumuData2019(const int nq=4, std::string beam="fhc", bool GetFromUPS=false)
TDirectory * dir
Definition: macro.C:5
const FitSinSqTheta23LowerOctant kFitSinSqTheta23LowerOctant
Definition: FitVars.cxx:17
void ResetToNominal()
Definition: SystShifts.cxx:141
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
virtual void SetTh23(const T &th23)=0
exit(0)
assert(nhit_max >=nhit_nbins)
Interface definition for fittable variables.
Definition: IFitVar.h:16
double GetPOT(bool isFHC=true)
std::vector< std::pair< const ISyst *, const ISyst * > > GetCorrelations(bool isNue, bool isFHC, bool ptExtrap)
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
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
Double_t signFunc(Double_t *x, Double_t *t)
And supporting functions:
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
std::pair< Spectrum *, double > cos
Float_t e
Definition: plot.C:35
#define for
Definition: msvc_pragmas.h:3
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
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)
void POTtag(bool both, bool FHCOnly, bool RHCOnly)
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
double GetX(int ndb=14, int db=1, int feb=0, int pix=0)
Definition: PlotOnMon.C:111
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
virtual void SetdCP(const T &dCP)=0
T asin(T number)
Definition: d0nt_math.hpp:60
std::pair< TH1D *, double > GetNueCosmics2019(std::string beam, bool GetFromUPS=false, bool NERSC=false)
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