reach_2018_dCPfractions.C
Go to the documentation of this file.
5 #include "CAFAna/Core/Progress.h"
6 //#include "CAFAna/Core/Utilities.h"
7 //#include "CAFAna/Experiment/CountingExperiment.h"
13 #include "CAFAna/Vars/FitVars.h"
15 
16 #include "OscLib/IOscCalc.h"
17 #include "Utilities/rootlogon.C"
18 
19 using namespace ana;
20 
21 #include "TCanvas.h"
22 #include "TFile.h"
23 #include "TGraph.h"
24 #include "TH2.h"
25 #include "TLatex.h"
26 #include "TLegend.h"
27 #include "TLine.h"
28 #include "TMath.h"
29 #include "TString.h"
30 
31 #include "TColor.h"
32 
33 #include <iostream>
34 #include <iomanip>
35 #include <vector>
36 
37 
38 ////////////////////////////////////////////////////////////////////////////////
40 {
41  std::vector <const IFitVar*> vars = {&kFitDeltaInPiUnits,
45  std::cout << "Oscillation parameters ";
46  for (const auto & v:vars){
47  std::cout << v->LatexName() << " " << v->GetValue(calc) << ", ";
48  }
50 }
51 
53  std::cout << "\n" << std::string(100, '*')<< std::endl;
54 }
55 
56 ////////////////////////////////////////////////////////////////////////////////
57 
58 struct SuperPred{
59  bool isNue;
61  double pot;
63 };
64 
65 std::vector <SuperPred> GetAllPredictions(const double potFHC,
66  const double potRHC,
67  const bool stats = false,
68  const bool nueOnly = false,
69  const bool numuOnly = false,
70  const std::string nueXP = "combo",
71  const std::string numuXP = "prop")
72 
73 {
74  std::vector <SuperPred> allpreds;
75  if(!numuOnly){
76  if(potFHC > 0) {
77  allpreds.push_back({true, "fhc", potFHC,
78  GetNuePrediction2018(nueXP, DefaultOscCalc(), !stats, "fhc",
79  true, false, true, false)});
80  }
81 
82  std::cout << "\n" << std::string(100, '*')<< std::endl;
83 
84  if(potRHC > 0) {
85  allpreds.push_back({true, "rhc", potRHC,
86  GetNuePrediction2018("prop", DefaultOscCalc(), !stats, "rhc",
87  true, false, true, false)});
88  }
89  } //end add nue predictions
90 
91  std::cout << "\n" << std::string(100, '*')<< std::endl;
92 
93  const int nnumu = 4;
94  if(!nueOnly){
95  if(potFHC > 0){
96  auto numu_preds = GetNumuPredictions2018(nnumu, !stats, "fhc");
97  for (int i = 0; i < nnumu; ++ i) {
98  allpreds.push_back({false,"fhc", potFHC, numu_preds[i]});
99  }
100  }
101 
102  std::cout << "\n" << std::string(100, '*')<< std::endl;
103 
104  if(potRHC > 0){
105  auto numu_preds = GetNumuPredictions2018(nnumu, !stats, "rhc");
106  for (int i = 0; i < nnumu; ++ i) {
107  allpreds.push_back({false,"rhc", potRHC, numu_preds[i]});
108  }
109  }
110  }//end add numu predictions
111 
112  return allpreds;
113 } //end get all predictions
114 
115 ////////////////////////////////////////////////////////////////////////////////
116 struct potCombo{
117  int year;
118  double potFHC;
119  double potRHC;
120 };
121 
122 std::vector <potCombo> GetPOTCombination(const TString potOption)
123 {
124  std::vector <potCombo> ret;
125  if(potOption == "default_short"){
126  ret = {
127  //{2017, 9.1, 0 },
128 // {2018, 9.1, 8.1},
129 // {2019, 11.6, 11.6},
130  {2019,12.61,12.61},
131 // {2020, 14.6, 14.6}
132  };
133  }
134 
135  if(potOption == "nominal_outdated"){
136  ret = { //{2016, 6.05,0.00},
137  {2017, 8.85,0.00},
138  {2018, 8.85,8.10},
139  {2019, 11.48,11.48},
140  {2020, 14.48,14.48},
141  {2021, 17.48,17.48},
142  {2022, 20.48,20.48},
143  {2023, 23.48,23.48},
144  {2024, 26.48,26.48}};
145  }
146  if(potOption == "beampower_updated"){
147  ret = { //{2016, 6.05,0.00},
148  {2017, 9.47,0.00},
149  {2018, 9.47,6.91},
150  {2019,12.61,12.61},
151  {2020,16.66,16.66},
152  {2021,21.22,21.22},
153  {2022,25.79,25.79},
154  {2023,30.85,30.85},
155  {2024,35.92,35.92} };
156  }
157 
158  std::cout << "\n POT combination is: " << potOption << "\n";
159  for(auto const & combo:ret){
160  std::cout << combo.year << ": \t"
161  << combo.potFHC<< " e20 POT FHC, \t"
162  << combo.potRHC<< " e20 POT RHC"
163  << std::endl;
164  }
165  return ret;
166 }
167 
168 TString GetPOTComboTitle(const TString potOption, bool stats = true){
169  TString ret;
170 
171  if(potOption=="nominal_outdated")
172  ret = "#splitline{DUMMY 2018 ana, stats only }{old nom POT combo } ";
173  if(potOption=="beampower_updated")
174  ret = "#splitline{2018 analysis techniques and}{projected beam exposure improvements} ";
175 
176  return ret;
177 }
178 
179 TLatex * PrintParamLabel (TString label){
180  TLatex * ltx = new TLatex();
181  ltx->SetTextAlign(13);
182  ltx->SetTextSize(0.04);
183  ltx->DrawLatexNDC(0.12,0.97,label.Data());
184  return ltx;
185 }
186 
187 ////////////////////////////////////////////////////////////////////////////////
188 TCanvas * PaintReachCanvas(TString param_label,
189  TString pot_option,
190  bool pot_label = true,
191  bool significance = false,
192  double maxY = 100,
193  TString yaxistitle = "")
194 {
195 
196  auto ret = new TCanvas;
197  gPad->SetTopMargin(0.10);
198 
199  auto combos = GetPOTCombination(pot_option);
200  int minX = combos.begin()->year;
201  int maxX = (combos.end()-1)->year;
202  int nBins = (combos.size() + 1) * 2;
203 
204  auto axes = new TH2D ("ax",";Year;Fraction of #delta_{CP} values (%)",
205  nBins, minX - 0.5, maxX + 0.5, 100, 0, maxY);
206  if(significance) axes->GetYaxis()->SetTitle("Significance #left(#sigma=#sqrt{#Delta#chi^{2}}#right)");
207  if(yaxistitle != "") axes->GetYaxis()->SetTitle(yaxistitle);
208 
209  axes->GetXaxis()->CenterTitle();
210  axes->GetYaxis()->CenterTitle();
211  axes->GetYaxis()->SetTitleOffset(0.6);
212  axes->Draw();
213 
214  PrintParamLabel(param_label);
215 
216  TLatex* ltx2 = new TLatex();
217  ltx2->SetTextSize(1.2*axes->GetXaxis()->GetLabelSize());
218  if(pot_label) ltx2->DrawLatexNDC(0.35, 0.16, GetPOTComboTitle(pot_option));
219 
220 // gPad->SetGridx();
221 // gPad->SetGridy();
222 
223  TLine * l1 = new TLine();
224  l1->SetLineStyle(3);
225  if(significance) l1->DrawLine(minX-0.5,3,maxX+0.5,3.); //3 sigma line
226 
227  SimulationSide();
228  return ret;
229 }
230 
231 double GetFractiondCP(TGraph * gr, const double level){
232 
233  //integration above level makes more sense. no time to learn
234  //transform the graph into a list of sigmas
235  std::vector <double>sigmas;
236  int nsteps = 1000;
237  for(int i=0;i<nsteps;++i){
238  sigmas.push_back(gr->Eval((double) i * 2./nsteps));
239  }
240  int numAbove = std::count_if(sigmas.begin(),sigmas.end(),
241  [level](double s){return s>=level;});
242 
243  return (double) 100* numAbove / sigmas.size();
244 }
245 
246 // New for 2018: Fitter hates us. Shiqi wrote this one; is it worth committing?
248 {
249 public:
250  virtual double GetValue(const osc::IOscCalcAdjustable* osc) const {
251  return util::sqr(sin(osc->GetTh23()));};
252  virtual void SetValue(osc::IOscCalcAdjustable* osc, double val) const {
253  osc->SetTh23(asin(sqrt(Clamp(val))));};
254  virtual std::string ShortName() const {return "ssth23MaxMix";}
255  virtual std::string LatexName() const {return "sin^{2}#theta_{23} MaxMix";}
256  virtual double LowLimit() const {return 0.49;}
257  virtual double HighLimit() const {return 0.51;}
258 };
260 
261 ////////////////////////////////////////////////////////////////////////////////
262 MultiExperiment * MakeFakeExperiment(const std::vector<SuperPred> allpreds,
263  osc::IOscCalcAdjustable* calcFake,
264  std::vector < std::pair <TH1D*, double> > cosmics = {},
265  std::vector <const IExperiment * > constraints = {},
266  bool verbose =false)
267 {
268 
269  bool hascosmics = cosmics.size()>0;
270  if(hascosmics && cosmics.size()!=allpreds.size()){
271  std::cerr << "\n Size of cosmics doesn't match predictions \n";
272  exit(1);
273  }
274  std::vector <const IExperiment * > expts;
275 
276  if(verbose) {
277  std::cout << "Creating Fake experiment \n";
278  PrintOscParameters(calcFake);
279  }
280 
281  for(int i = 0; i < (int) allpreds.size(); ++i){
282 
283  auto thisdata = GetFakeData (allpreds[i].prediction, calcFake, allpreds[i].pot,
284  hascosmics?cosmics[i].first:0);
285  if(allpreds[i].pot > 0){
286  if(hascosmics){
287  expts.push_back(new SingleSampleExperiment(allpreds[i].prediction, *thisdata,
288  cosmics[i].first, cosmics[i].second) );
289  }
290  else {
291  expts.push_back(new SingleSampleExperiment(allpreds[i].prediction, *thisdata) ) ;
292  }
293 
294  if(verbose){
295  auto temp = thisdata->ToTH1(thisdata->POT());
296  std::cout << "Added experiment "
297  << allpreds[i].pot << " POT "
298  << allpreds[i].flux << " "
299  << (allpreds[i].isNue?"nue ":"numu ")
300  << temp->Integral()
301  << " events\n";
302  delete temp;
303  }
304  }
305  }
306 
307  // Add additional constraints and create expt
308  if(verbose) {
309  std::cout << "Returning Multiexperiment from " << expts.size()
310  << " SingleSample experiments and "
311  << constraints.size() << " external constraints" << std::endl;
312  }
313  expts.insert(expts.end(), constraints.begin(), constraints.end());
314 
315  MultiExperiment * expt = new MultiExperiment (expts);
316 
317  //Maybe set sys correlations here? or add a new function
318  return expt;
319 }
320 
321 
322 ////////////////////////////////////////////////////////////////////////////////
323 void PlotDeltaFractionVsYear(TString meas,
324  std::vector <TString> hierarchies,
325  std::vector <TString> thetas,
326  std::vector <TString> combo_names,
327  std::vector <Int_t > colors,
328  bool react = true,
329  std::map < TString, TString > ssth23_label = {},
330  std::map < TString, TString> combo_label={})
331 {
332  double limit = 3.84; //dchisq for 95%o
333 
334  TLegend * leg = new TLegend (0.12, 0.6, 0.5, 0.89);
335  leg->SetFillStyle(0);
336  leg->SetMargin(1.3*leg->GetMargin());
337  if( combo_names.size()>1 ) leg->SetNColumns(2);
338  leg->SetTextSize(0.04);
339 
340  TString header = "";//"#splitline{NOvA joint #nu_{e}+#nu_{#mu}}{";
341  if(meas == "hie") header += "Hierarchy at 95% CL";
342  if(meas == "cpv") header += "CPV at 95% CL";
343 // header+="}";
344  leg->SetHeader(header);
345 
346 
347  for(uint i = 0; i < combo_names.size(); ++i){
348 
349  auto combos = GetPOTCombination (combo_names[i]);
350 
351  int col = 0;
352  for (auto const & hie: hierarchies){
353  for (auto const & theta:thetas){
354 
355  TString fname = ("dCPfraction/hist_" + meas+ "_" + hie + "_" +
356  theta + "_" + (react ? "react_" : "") +
357  combo_names[i] + ".root") ;
358 
359  auto file = new TFile (fname,"read");
360  if(!file) {std::cerr << "No file " << fname << std::endl; exit(1);}
361  auto gr_frac = new TGraph();
362 
363  for(auto & combo:combos){
364 
365  auto gr = (TGraph*) file->Get(meas + "_" + hie + "_" + theta + "_" +
366  (react ? "react_" : "") +
367  combo_names[i] +
368  "_" + std::to_string(combo.year) + "/gr_dchisq");
369 
370  if(!gr) {std::cerr << "Problem w graph \n"; continue;}
371  auto thisfrac = GetFractiondCP(gr, limit);
372  gr_frac->SetPoint(gr_frac->GetN(), combo.year, thisfrac);
373  }
374  gr_frac->SetLineWidth(2);
375  gr_frac->SetLineStyle(i+1);
376  gr_frac->Draw("L");
377  gr_frac->SetLineColor(colors[col]);
378  ++col;
379  leg->AddEntry(gr_frac, hie + " " + ssth23_label[theta] + " " +
380  combo_label[combo_names[i]],"l");
381  }//end ssth23
382  }//end hie
383  }//end combo names
384  leg->Draw();
385  gPad->RedrawAxis();
386 }
387 
388 
390  std::vector <TString> hierarchies,
391  std::vector <TString> thetas,
392  std::vector <TString> combo_names,
393  std::vector <Int_t > colors,
394  bool react = true,
395  std::map < TString, TString > ssth23_label = {},
396  std::map < TString, TString> combo_label={})
397 {
398  gPad->SetBottomMargin(0.12);
399 
400  TString head= "";
401  double maxx = 6.;
402  if(meas == "hie") head = "Hierarchy resolution";
403  if(meas == "cpv") { head = "CP violation"; maxx =3.5;}
404  auto leg2 = new TLegend(0.6,0.6,0.9,0.9);
405  leg2->SetFillStyle(0);
406  leg2->SetHeader(head);
407 
408  auto ax2 = new TH2D("ax2",
409  ";Significance #left(#sigma=#sqrt{#Delta#chi^{2}}#right);Fraction of #delta_{CP} values (%)",
410  100,0,maxx,100,0,100);
411  ax2->Draw();
412  ax2->GetXaxis()->CenterTitle();
413  ax2->GetYaxis()->CenterTitle();
414  ax2->GetYaxis()->SetTitleOffset(0.8);
415 
416  int col =0;
417  for (auto const & hie: hierarchies){
418  for (int th = 0; th < 2; ++th){
419 
420  TString fname = "dCPfraction/hist_" + meas + "_" + hie +"_"+thetas[th] + "_" +
421  (react? "react_":"") + combo_names[0] + ".root";
422  auto file = new TFile (fname,"read");
423  if(!fname) exit(1);
424 
425  for(int year:{2024}){
426  auto gr = (TGraph*) file->Get(meas + "_"+ hie +"_" + thetas[th] + "_" +
427  (react?"react_":"") + combo_names[0] +
428  "_" + std::to_string(year) + "/gr_dchisq");
429  if(!gr) {std::cerr << "Problem w graph \n"; continue;}
430 
431  auto gr_frac = new TGraph();
432 
433  for(double lim = 0; lim < 6; lim += 0.05 ){
434  auto thisfrac = GetFractiondCP(gr, lim*lim);
435  gr_frac->SetPoint(gr_frac->GetN(), lim,thisfrac);
436  }
437  gr_frac->SetLineWidth(2);
438  gr_frac->SetLineStyle(th+1);
439  gr_frac->Draw("L");
440  gr_frac->SetLineColor(colors[col]);
441  leg2->AddEntry (gr_frac, hie + " " + ssth23_label[thetas[th]], "l");
442 // if (th==0) leg2->AddEntry(gr_frac, std::to_string(year).c_str(), "l");
443 // ++col;
444  }//end years
445  }//end th
446  ++col;
447  }//end hie
448 /* auto dummy = new TGraph();
449  dummy->SetLineColor(kGray+1);
450  dummy->SetLineWidth(2);
451  leg2->AddEntry(dummy->Clone(),"sin^{2}#theta_{23}=0.6","l");
452  dummy->SetLineStyle(2);
453  leg2->AddEntry(dummy->Clone(),"sin^{2}#theta_{23}=0.4","l");
454 */
455  leg2->SetTextSize(ax2->GetXaxis()->GetLabelSize());
456  leg2->Draw();
457 // gPad->SetGrid();
458  gPad->RedrawAxis();
459 }
460 ////////////////////////////////////////////////////////////////////////////////
462  std::vector <TString> hierarchies,
463  std::vector <double> deltas,
464  std::vector <TString> thetas,
465  std::vector <TString> combo_names,
466  std::vector <Int_t > colors,
467  bool react = true,
468  std::map < double, TString> delta_label = {}
469  )
470 {
471  TString head= "";
472  if(meas == "hie") head = "Hierarchy resolution";
473  if(meas == "cpv") head = "CP violation";
474  if(meas == "oct") head = "Octant determination";
475  if(meas == "mix") head = "Max. mixing rejection";
476 
477  auto leg3 = new TLegend(0.1,0.6,0.5,0.89);
478  leg3->SetHeader(head);
479  leg3->SetTextSize(0.04);
480 
481  for(uint i = 0; i < combo_names.size(); ++i){
482 
483  auto combos = GetPOTCombination (combo_names[i]);
484  int col = 0;
485  for (auto const & hie: hierarchies){
486  for(double delta:deltas){
487 
488  auto gr_max = new TGraph();
489  auto gr_min = new TGraph();
490 
491  for(auto & combo:combos){
492 
493  double maxchi = -999.; double minchi = 999.;
494 
495  for (auto const & theta:thetas){
496  TString fname = "dCPfraction/hist_" + meas + "_" + hie + "_" + theta + "_" +
497  (react ? "react_":"") + combo_names[i] + ".root";
498 
499  TFile file (fname,"read");
500  if(file.IsZombie()) {std::cerr << "No file " << fname << std::endl; exit(1);}
501 
502  auto gr = (TGraph*) file.Get(meas + "_" + hie + "_" + theta + "_" +
503  (react ? "react_":"") + combo_names[i] +
504  "_" + std::to_string(combo.year) + "/gr_dchisq");
505 
506  if(!gr) {std::cerr << "Problem w graph \n"; continue;}
507 
508  maxchi = std::max(maxchi, gr->Eval(delta));
509  minchi = std::min(minchi, gr->Eval(delta));
510  } //end thetas
511  gr_max->SetPoint(gr_max->GetN(),combo.year, sqrt(maxchi));
512  gr_min->SetPoint(gr_min->GetN(),combo.year, sqrt(minchi));
513  }//end years
514  auto tempN = gr_min->GetN();
515  for(int i = 0; i < tempN; ++i ){
516  gr_max->SetPoint(gr_max->GetN(),
517  gr_min->GetX()[tempN - i - 1],
518  gr_min->GetY()[tempN - i - 1]);
519  }
520 
521  gr_max->Draw("f");
522  gr_max->SetFillColorAlpha(colors[col],.33);
523  ++col;
524  gr_max->SetLineColor(gr_max->GetFillColor());
525  if(delta_label.size()==0) {
526  leg3->AddEntry(gr_max, hie + TString::Format(" #delta_{CP}= %.2f#pi",delta),"f");}
527  else{
528  leg3->AddEntry(gr_max, hie + " #delta_{CP}="+delta_label[delta],"f");
529  }
530  }//end deltas
531  }//end hie
532  }//end combo names
533  leg3->SetFillStyle(0);
534  leg3->Draw();
535 }
536 
537 
538 
539 ////////////////////////////////////////////////////////////////////////////////
541  std::vector <TString> hierarchies,
542  std::vector <TString> thetas,
543  std::vector <TString> combo_names,
544  std::vector <Int_t > colors,
545  bool react = true,
546  std::map < TString, TString > ssth23_label = {})
547 {
548  TString head= "";
549  if(meas == "hie") head = "Hierarchy resolution";
550  if(meas == "cpv") head = "CP violation";
551  if(meas == "oct") head = "Octant determination";
552  if(meas == "mix") head = "Max. mixing rejection";
553 
554  auto leg3 = new TLegend(0.1,0.6,0.5,0.88);
555  leg3->SetHeader(head);
556  leg3->SetTextSize(0.04);
557 
558  for(uint i = 0; i < combo_names.size(); ++i){
559 
560  auto combos = GetPOTCombination (combo_names[i]);
561  int col = 0;
562  for (auto const & hie: hierarchies){
563  for (auto const & theta:thetas){
564 
565  TString fname = "dCPfraction/hist_" + meas + "_" + hie + "_" + theta + "_" +
566  (react ? "react_":"") + combo_names[i] + ".root";
567 
568  TFile file (fname,"read");
569  if(file.IsZombie()) {std::cerr << "No file " << fname << std::endl; exit(1);}
570 
571  auto gr_max = new TGraph();
572  auto gr_min = new TGraph();
573 
574  for(auto & combo:combos){
575 
576  double maxchi = -999.; double minchi = 999.;
577 
578  auto gr = (TGraph*) file.Get(meas + "_" + hie +"_" + theta + "_" +
579  (react ? "react_":"") + combo_names[i] +
580  "_" + std::to_string(combo.year) + "/gr_dchisq");
581 
582  if(!gr) {std::cerr << "Problem w graph \n"; continue;}
583 
584  minchi = TMath::MinElement(gr->GetN(),gr->GetY());
585  maxchi = TMath::MaxElement(gr->GetN(),gr->GetY());
586 
587  gr_max->SetPoint(gr_max->GetN(),combo.year, sqrt(maxchi));
588  gr_min->SetPoint(gr_min->GetN(),combo.year, sqrt(minchi));
589  } //end years
590  auto tempN = gr_min->GetN();
591  for(int i = 0; i < tempN; ++i ){
592  gr_max->SetPoint(gr_max->GetN(),
593  gr_min->GetX()[tempN - i - 1],
594  gr_min->GetY()[tempN - i - 1]);
595  }
596 
597  gr_max->Draw("f");
598  gr_max->SetFillColorAlpha(colors[col],.33);
599  ++col;
600  gr_max->SetLineColor(gr_max->GetFillColor());
601  if(ssth23_label.size()==0) {
602  leg3->AddEntry(gr_max, hie + " " + theta, "f");}
603  else{
604  leg3->AddEntry(gr_max, hie + " "+ssth23_label[theta],"f");
605  }
606  }//end theta
607  }//end hier
608  }//end combo names
609  leg3->SetFillStyle(0);
610  leg3->Draw();
611 }
612 
613 
614 ////////////////////////////////////////////////////////////////////////////////
615 
616 void reach_2018_dCPfractions(const bool createFile = false,
617  const TString options = "",
618  const TString pot_combo = "default_short")
619 {
620  // Shared info
621  TString filename = "dCPfraction/hist_"+ options + "_" + pot_combo + ".root";
622  PrintSeparator();
623 
624  auto combos = GetPOTCombination(pot_combo);
625  PrintSeparator();
626 
627  //////////////////////////////////////////////////////////////////////
628  bool isNH = options.Contains("NH");
629  bool isIH = options.Contains("IH");
630  bool isUO = options.Contains("UO");
631  bool isLO = options.Contains("LO");
632  bool isMAX = options.Contains("MAX");
633 
634  assert(isNH || isIH || !createFile);
635  assert(isUO || isLO || isMAX || !createFile);
636 
637  double ssth23 = (isUO ? 0.6 : 0.4);
638  if(isMAX) ssth23 = 0.5;
639 
640  double dmsq32 = 2.50e-3;
641  if(!isNH) dmsq32 = -2.55e-3;
642 
643  double dcp_pi = 0.2;
644 
645  bool run_hie = options.Contains("hie");
646  bool run_cpv = options.Contains("cpv");
647  bool run_oct = options.Contains("oct");
648  bool run_mix = options.Contains("mix");
649 
650  assert(run_hie || run_cpv || run_oct || run_mix);
651 
652  bool reactor = options.Contains("react");
653 
654  if( (run_mix || run_oct) && isMAX ) {
655  std::cerr << "No wrong octant or rejection of maximal mixing if MAX is true"
656  << std::endl;
657  exit(1);
658  }
659 
660  //////////////////////////////////////////////////////////////////////
661 
662  if(createFile){
663 
664  auto file = new TFile(filename, "recreate");
665  auto preds = GetAllPredictions(1E20, 1E20, true);
666  PrintSeparator();
667  //////////////////////////////////////////////////////////////////////
668  //Setup for hierarchy,cpv, dcp fractions
669  double dstep = 2./32.;
670 
671  std::vector <const IFitVar *> fitvars_hie = {&kFitDeltaInPiUnits,
673  &kFitDmSq32
674  };
675  std::vector <const IFitVar *> fitvars_cpv = {
677  &kFitDmSq32
678  };
679  std::vector <const IFitVar *> fitvars_oct = {&kFitDeltaInPiUnits,
680  isUO ?
683  &kFitDmSq32
684  };
685  std::vector <const IFitVar *> fitvars_mix = {&kFitDeltaInPiUnits,
687  &kFitDmSq32
688  };
689 
690  std::map <const IFitVar*, std::vector <double> > seeds_hie =
691  {
692  {&kFitDeltaInPiUnits, {0.,0.5,1.,1.5}},
693  {&kFitSinSqTheta23, {0.3,0.7}}
694  };
695  std::map <const IFitVar*, std::vector <double> > seeds_cpv =
696  {
697  {&kFitSinSqTheta23, {0.3,0.7}}
698  };
699  std::map <const IFitVar*, std::vector <double> > seeds_oct =
700  {
701  {&kFitDeltaInPiUnits, {0.,0.5,1.,1.5}}
702  };
703  std::map <const IFitVar*, std::vector <double> > seeds_mix =
704  {
705  {&kFitDeltaInPiUnits, {0.,0.5,1.,1.5}},
706  };
707  if(isUO) seeds_mix.insert({kFitSinSqTheta23LowerOctant,{0.49}});
708  if(isLO) seeds_mix.insert({kFitSinSqTheta23UpperOctant,{0.51,0.6}});
709 
710  std::vector <const IExperiment* > constraints = {};
711  if(reactor) {
712  for (auto fv:{&fitvars_hie, &fitvars_cpv, &fitvars_oct, &fitvars_mix}){
713  fv->push_back(&kFitSinSq2Theta13);
714  }
715  constraints.push_back(WorldReactorConstraint2017());
716  }
717  else {
718  std::cerr << "\nWARNING: Not using reactor constraint\n";
719  }
720  //////////////////////////////////////////////////////////////////////
721  for (auto const & combo:combos){
722  PrintSeparator();
723  std::cout << "\n\n Start POT combination "
724  << combo.year << ": \t"
725  << combo.potFHC<< " e20 POT FHC, \t"
726  << combo.potRHC<< " e20 POT RHC \n"
727  << std::endl;
728 
729  //Change the pot in the super predictions
730  for(auto & pred:preds){
731  if(pred.flux=="fhc") pred.pot = combo.potFHC * 1E20;
732  if(pred.flux=="rhc") pred.pot = combo.potRHC * 1E20;
733  }
734  // Create an identifier and folder based on the pot combo
735  // Store the POT info just in case
736  TString dirYear = ( options + "_" + pot_combo + "_" +
737  std::to_string(combo.year) );
738 
739  auto dir = file->mkdir(dirYear);
740  dir->WriteObject(&combo,"potCombo");
741 
742  // Store the fake data calculator just in case - dCP will change
743  auto calc_fake = DefaultOscCalc();
744  SetFakeCalc(calc_fake, ssth23, dmsq32, dcp_pi);
745  ana::SaveTo(* (osc::IOscCalc*) calc_fake, dir, "calc");
746 
747  //One graph = one dCP loop. Set steps earlier
748  TGraph* gr = new TGraph();
749 
750  Progress prog ("Starting dCP loop");
751 
752  for(double del = 0; del < 2; del += dstep){
753  prog.SetProgress(del/2);
754 
755  SetFakeCalc(calc_fake, ssth23, dmsq32, del);
756 
757  // Create the experiment
758  std::vector < std::pair <TH1D*, double > > cosmics = {};
759 
760  //Some verbosity to check the fake experiments, but not all
761  bool verb = (del == 0);
762  auto thisexpt = MakeFakeExperiment(preds, calc_fake,
763  cosmics, constraints, verb);
764 
765  //maybe add correlations here
766 
767  //placeholder
768  std::vector <const ISyst* > mySysts = {};
769 
770  double minchi_true = 0;
771  double minchi_wrong = 1e20;
772  SystShifts systSeed = SystShifts::Nominal();
773  auto calc = DefaultOscCalc();
774 
775  if(run_hie){
776  // Get fits in wrong hierarchy
777  MinuitFitter fit_hie(thisexpt, fitvars_hie, mySysts);
778  calc->SetDmsq32(-1. * calc_fake->GetDmsq32()); //todo
779  minchi_wrong = fit_hie.Fit(calc,systSeed, seeds_hie, {},
780  verb ? IFitter::kVerbose : IFitter::kQuiet)->EvalMetricVal();
781  }
782  if(run_cpv){
783  MinuitFitter fit_cpv(thisexpt, fitvars_cpv, mySysts);
784  for(int hie:{-1,1}) { //test both hierarchies
785  for(int dcp:{0,1}) { //test cp conserving values
786  systSeed.ResetToNominal();
787  calc = DefaultOscCalc();
788  calc->SetDmsq32(hie * calc_fake->GetDmsq32());
789  calc->SetdCP(dcp*M_PI);
790  minchi_wrong = std::min(minchi_wrong,
791  fit_cpv.Fit(calc,systSeed, seeds_cpv, {},
792  verb ? IFitter::kVerbose :
793  IFitter::kQuiet)->EvalMetricVal());
794  }//end dcp
795  }//end hie
796  } //end cpv
797 
798  if(run_oct){
799  MinuitFitter fit_oct(thisexpt, fitvars_oct, mySysts);
800  for(int hie:{-1,1}) { //test both hierarchies
801  systSeed.ResetToNominal();
802  calc = DefaultOscCalc();
803  calc->SetDmsq32(hie * calc_fake->GetDmsq32());
804  calc->SetTh23(isUO? asin(sqrt(0.4)): asin(sqrt(0.6)));
805  //wrong octant is set through seeds and fitvars
806  minchi_wrong = std::min(minchi_wrong,
807  fit_oct.Fit(calc,systSeed, seeds_oct, {},
808  verb ? IFitter::kVerbose :
809  IFitter::kQuiet)->EvalMetricVal());
810  }//end hie
811  } //end oct
812 
813  if(run_mix){
814  MinuitFitter fit_mix(thisexpt, fitvars_mix, mySysts);
815  for(int hie:{-1,1}) { //test both hierarchies
816  systSeed.ResetToNominal();
817  calc = DefaultOscCalc();
818  calc->SetDmsq32(hie * calc_fake->GetDmsq32());
819  calc->SetTh23(asin(sqrt(0.5)));
820  minchi_wrong = std::min(minchi_wrong,
821  fit_mix.Fit(calc,systSeed, seeds_mix, {},
822  verb ? IFitter::kVerbose :
823  IFitter::kQuiet)->EvalMetricVal());
824  }//end hie
825  } //end mix
826 
827  // Check best fit overall sometimes; todo: some cases actually want
828  // to store, not set to 0
829  if(verb && 1){
830  MinuitFitter fit_best(thisexpt, fitvars_hie, mySysts);
831  for(int hie:{-1,1}) {
832  systSeed.ResetToNominal();
833  calc = DefaultOscCalc();
834  calc->SetDmsq32(hie * calc_fake->GetDmsq32());
835  minchi_true = fit_best.Fit(calc,systSeed, seeds_hie)->EvalMetricVal();
836  std::cout << "\nFound true for " << (hie > 0 ? "NH" : "IH")
837  << " dCP/pi" << del << " minchi " << minchi_true
838  << std::endl;
839  minchi_true = 0; //bad bad hack hack
840  }
841  }
842  gr->SetPoint(gr->GetN(),del,minchi_wrong - minchi_true);
843  }//end dCP loop
844  prog.Done();
845  dir->cd();
846  gr->Write("gr_dchisq");
847  }//end combos
848  std::cout << "Wrote graphs to file " << filename << std::endl;
849  }//end create file
850 
851 ////////////////////////////////////////////////////////////////////////////////
852  else{ //make plots
853 
854  std::vector <TString> thetas = {};
855  if(isUO || 1) thetas.push_back("UO");
856  if(isMAX || 1) thetas.push_back("MAX");
857  if(isLO || 1) thetas.push_back("LO");
858 
859  std::map < TString, TString > ssth23_label = {{"UO","sin^{2}#theta_{23}=0.6"},
860  {"MAX","sin^{2}#theta_{23}=0.5"},
861  {"LO","sin^{2}#theta_{23}=0.4"}
862  };
863 
864  std::vector <TString> combo_names = {}; //sue me
865  if(0) combo_names.push_back("nominal_outdated");
866  if(1) combo_names.push_back("beampower_updated");
867 
868  std::map < TString, TString > combo_label = {{"beampower_updated",""}};
869 
870  std::vector <Int_t > Colors = {
871  //violet blue cyan green yellow orange red magenta
872  TColor::GetColor("#6c71c4"), TColor::GetColor("#268bd2"),
873  TColor::GetColor("#2aa198"), TColor::GetColor("#859900"),
874  TColor::GetColor("#b58900"), TColor::GetColor("#cb4b16"),
875  TColor::GetColor("#dc322f"), TColor::GetColor("#d33682")
876  };
877  std::vector <Int_t > Colors_alt = {
878  //cyan violet yellow magenta blue red orange green
879  TColor::GetColor("#2aa198"),
880  TColor::GetColor("#6c71c4"),
881  TColor::GetColor("#b58900"),
882  TColor::GetColor("#d33682"),
883  TColor::GetColor("#268bd2"),
884  TColor::GetColor("#dc322f"),
885  TColor::GetColor("#cb4b16"),
886  TColor::GetColor("#859900"),
887  };
888 
889 
890  ////////////////////////////////////////////////////////////////////////////////
891  // dCP fraction plots
892  ////////////////////////////////////////////////////////////////////////////////
893  TString labelShort = TString::Format("|#Deltam^{2}_{32}|=2.5#times10^{-3}eV^{2}, sin^{2}2#theta_{13}=0.082"); //todo: not generic
894  //HIE
895  //Having these here gives me more control. Should add a check that canvas exist or sth
896  PaintReachCanvas(labelShort,combo_names[0],false,false,100);
897  PlotDeltaFractionVsYear("hie", {"NH"}, thetas, combo_names,
898  Colors_alt,
899  reactor, ssth23_label, combo_label);
900  gPad->Print("reach2018_hie_dCPfrac_year_trueNH.pdf");
901 
902  PaintReachCanvas(labelShort,combo_names[0],false,false,100);
903  PlotDeltaFractionVsYear("hie", {"IH"}, thetas, combo_names,
904  Colors_alt,
905  reactor, ssth23_label, combo_label);
906  gPad->Print("reach2018_hie_dCPfrac_year_trueIH.pdf");
907 
908  // CPV
909  PaintReachCanvas(labelShort,combo_names[0],false,false,60);
910 
911  PlotDeltaFractionVsYear("cpv", {"NH"}, thetas, combo_names,
912  Colors_alt, reactor, ssth23_label, combo_label);
913  gPad->Print("reach2018_cpv_dCPfrac_year_trueNH.pdf");
914 
915  PaintReachCanvas(labelShort,combo_names[0],false,false,60);
916  PlotDeltaFractionVsYear("cpv", {"IH"}, thetas, combo_names,
917  Colors_alt, reactor, ssth23_label, combo_label);
918  gPad->Print("reach2018_cpv_dCPfrac_year_trueIH.pdf");
919 
920 
921  ////////////////////////////////////////////////////////////////////////////////
922  // Delta fraction vs significance for year 2024
923  TString label2024 = TString::Format("#Deltam^{2}_{32}=2.5#times10^{-3}eV^{2}, sin^{2}2#theta_{13}=0.082. #kern[1.5]{ } 36#times10^{20} POT(#nu)+ 36#times10^{20} POT(#bar{#nu})");
924 
925  TLine * l1 = new TLine();
926  l1->SetLineStyle(3);
927 
928  //HIE
929  new TCanvas();
930  PlotDeltaFractionVsSignificance("hie", {"NH","IH"}, {"UO","LO"}, combo_names,
931  {Colors[1],Colors[6]},
932  reactor, ssth23_label, combo_label);
933  PrintParamLabel(label2024);
934  SimulationSide();
935  l1->DrawLine(3,0,3,100); //3 sigma line
936  gPad->Print("reach2018_hie_frac_sigma.pdf");
937 
938  //CPV
939  new TCanvas();
940  PlotDeltaFractionVsSignificance("cpv", {"NH","IH"}, {"UO","LO"}, combo_names,
941  {Colors[1],Colors[6]},
942  reactor, ssth23_label, combo_label);
943  PrintParamLabel(label2024);
944  SimulationSide();
945  l1->DrawLine(3,0,3,60); //3 sigma line
946  gPad->Print("reach2018_cpv_frac_sigma.pdf");
947 
948  //////////////////////////////////////////////////////////////////////////
949  // XX vs year plots
950  // Fixed delta: hie and cpv
951  std::vector <double> deltas = {1.5,1.,0.,0.5};
952  std::map <double, TString> delta_label = {{1.5,"3#pi/2"},{0.5," #pi/2"},
953  {1.,"#pi"},{0.,"0"} };
954 
955  TString toplabel = TString::Format("sin^{2}#theta_{23} #in #left{0.4, 0.5, 0.6#right}, |#Deltam^{2}_{32}|=2.5#times10^{-3}eV^{2}, sin^{2}2#theta_{13}=0.082");
956 
957  //HIE
958  PaintReachCanvas(toplabel,combo_names[0],true,true,5);
959  PlotSignificanceVsYear_FixDelta("hie", {"NH"}, deltas, thetas, combo_names,
960  Colors_alt, reactor, delta_label);
961  gPad->Print("reach2018_hie_sigma_year_trueNH.pdf");
962 
963  PaintReachCanvas(toplabel,combo_names[0],true,true,5);
964  PlotSignificanceVsYear_FixDelta("hie", {"IH"}, deltas, thetas, combo_names,
965  Colors_alt, reactor, delta_label);
966  gPad->Print("reach2018_hie_sigma_year_trueIH.pdf");
967 
968  //CPV
969  PaintReachCanvas(toplabel,combo_names[0],true,true,5);
970  PlotSignificanceVsYear_FixDelta("cpv", {"NH"}, {1.5,0.5}, thetas, combo_names,
971  Colors_alt, reactor, delta_label);
972  SimulationSide();
973  gPad->Print("reach2018_cpv_sigma_year_trueNH.pdf");
974 
975  PaintReachCanvas(toplabel,combo_names[0],true,true,5);
976  PlotSignificanceVsYear_FixDelta("cpv", {"IH"}, {1.5,0.5}, thetas, combo_names,
977  Colors_alt, reactor, delta_label);
978  SimulationSide();
979  gPad->Print("reach2018_cpv_sigma_year_trueIH.pdf");
980 
981  // Fixed ssth23: oct and maxmix
982  TString toplabel2 = TString::Format("#delta_{CP} #in #left[0,2#pi#right), |#Deltam^{2}_{32}=2.5#times10^{-3}eV^{2}|, sin^{2}2#theta_{13}=0.082");
983 
984  //OCT
985  PaintReachCanvas(toplabel2,combo_names[0],true,true,5);
986  PlotSignificanceVsYear_FixTh23("oct", {"NH","IH"}, {"UO","LO"}, combo_names,
987  Colors_alt, reactor, ssth23_label);
988 
989  gPad->Print("reach2018_oct_sigma_year.pdf");
990 
991  //MIX
992  PaintReachCanvas(toplabel2,combo_names[0],true,true,7);
993  PlotSignificanceVsYear_FixTh23("mix", {"NH","IH"}, {"UO","LO"}, combo_names,
994  Colors_alt, reactor, ssth23_label);
995  SimulationSide();
996  gPad->Print("reach2018_mix_sigma_year.pdf");
997  }//end make plots
998 } //end all
Spectrum * GetFakeData(const IPrediction *pred, osc::IOscCalc *calc, const double pot, TH1D *cosmics=0)
std::vector< potCombo > GetPOTCombination(const TString potOption)
T max(const caf::Proxy< T > &a, T b)
int nBins
Definition: plotROC.py:16
void SimulationSide()
Definition: rootlogon.C:137
void PrintOscParameters(const osc::IOscCalcAdjustable *calc)
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const IPrediction * prediction
void reach_2018_dCPfractions(const bool createFile=false, const TString options="", const TString pot_combo="default_short")
double ssth23
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
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
TLatex * PrintParamLabel(TString label)
const IConstrainedFitVar * kFitSinSqTheta23MaxMix
const int nsteps
T sqrt(T number)
Definition: d0nt_math.hpp:156
const FitSinSqTheta23UpperOctant kFitSinSqTheta23UpperOctant
Definition: FitVars.cxx:16
General interface to oscillation calculators.
Definition: StanTypedefs.h:22
void PlotDeltaFractionVsYear(TString meas, std::vector< TString > hierarchies, std::vector< TString > thetas, std::vector< TString > combo_names, std::vector< Int_t > colors, bool react=true, std::map< TString, TString > ssth23_label={}, std::map< TString, TString > combo_label={})
TString GetPOTComboTitle(const TString potOption, bool stats=true)
int stats(TString inputFilePath, Int_t firstRun, Int_t lastRun, Float_t thresh, TString myDet)
Definition: stats.C:13
OStream cerr
Definition: OStream.cxx:7
virtual double LowLimit() const
virtual T GetTh23() const
Definition: IOscCalc.h:94
double maxY
Definition: plot_hist.C:10
string filename
Definition: shutoffs.py:106
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)
static SystShifts Nominal()
Definition: SystShifts.h:34
constraints
Definition: mkDefs.py:147
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
void PrintSeparator()
#define M_PI
Definition: SbMath.h:34
int Clamp(int x)
Eigen::Matrix< T, Eigen::Dynamic, 1 > head(const Eigen::Matrix< T, Eigen::Dynamic, 1 > &v, size_t n)
Definition: head.hpp:24
void PlotSignificanceVsYear_FixTh23(TString meas, std::vector< TString > hierarchies, std::vector< TString > thetas, std::vector< TString > combo_names, std::vector< Int_t > colors, bool react=true, std::map< TString, TString > ssth23_label={})
const char * label
std::vector< SuperPred > GetAllPredictions(const double potFHC, const double potRHC, const bool stats=false, const bool nueOnly=false, const bool numuOnly=false, const std::string nueXP="combo", const std::string numuXP="prop")
const XML_Char * s
Definition: expat.h:262
virtual double GetValue(const osc::IOscCalcAdjustable *osc) const
expt
Definition: demo5.py:34
int colors[6]
Definition: tools.h:1
Int_t col[ntarg]
Definition: Style.C:29
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
#define pot
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
const std::map< std::pair< std::string, std::string >, Variable > vars
void PlotDeltaFractionVsSignificance(TString meas, std::vector< TString > hierarchies, std::vector< TString > thetas, std::vector< TString > combo_names, std::vector< Int_t > colors, bool react=true, std::map< TString, TString > ssth23_label={}, std::map< TString, TString > combo_label={})
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
void SetFakeCalc(osc::IOscCalcAdjustable *calc, double ssth23=-5, double dmsq32=-5, double dCP_Pi=-5)
virtual std::string LatexName() const
Oscillation probability calculators.
Definition: Calcs.h:5
MultiExperiment * MakeFakeExperiment(const std::vector< SuperPred > allpreds, osc::IOscCalcAdjustable *calcFake, std::vector< std::pair< TH1D *, double > > cosmics={}, std::vector< const IExperiment * > constraints={}, bool verbose=false)
virtual std::string ShortName() const
OStream cout
Definition: OStream.cxx:6
Combine multiple component experiments.
T sin(T number)
Definition: d0nt_math.hpp:132
void SetProgress(double frac)
Update the progress fraction between zero and one.
Definition: Progress.cxx:41
double dmsq32
TDirectory * dir
Definition: macro.C:5
const FitSinSqTheta23LowerOctant kFitSinSqTheta23LowerOctant
Definition: FitVars.cxx:17
void ResetToNominal()
Definition: SystShifts.cxx:143
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
virtual void SetTh23(const T &th23)=0
exit(0)
virtual double HighLimit() const
void PlotSignificanceVsYear_FixDelta(TString meas, std::vector< TString > hierarchies, std::vector< double > deltas, std::vector< TString > thetas, std::vector< TString > combo_names, std::vector< Int_t > colors, bool react=true, std::map< double, TString > delta_label={})
TFile * file
Definition: cellShifts.C:17
assert(nhit_max >=nhit_nbins)
A simple ascii-art progress bar.
Definition: Progress.h:9
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 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)
TCanvas * PaintReachCanvas(TString param_label, TString pot_option, bool pot_label=true, bool significance=false, double maxY=100, TString yaxistitle="")
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
T min(const caf::Proxy< T > &a, T b)
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
void Done()
Call this when action is completed.
Definition: Progress.cxx:90
virtual void SetValue(osc::IOscCalcAdjustable *osc, double val) const
double GetFractiondCP(TGraph *gr, const double level)
T asin(T number)
Definition: d0nt_math.hpp:60
Compare a single data spectrum to the MC + cosmics expectation.
static constexpr Double_t year
Definition: Munits.h:185
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17
enum BeamMode string
unsigned int uint