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